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ABSTRACT 

Passive  microwave  satellite  data  provides  all-weather,  continuous 
observations  of  the  polar  ice  cover.  Multichannel  radiometer  data  sets  extending 
back  to  1978  are  readily  available  and  distributed  on  CD-ROM  media.  Current 
algorithms  can  extract  information  on  ice  cover,  concentration  and  type.  Two 
applications  of  these  data  sets  for  polar  studies  are  explored:  1)  the  possible 
detection  of  altered  ice  morphology  due  to  anomalous  tidal  currents  in  the 
Yermak  Plateau  region;  and  2)  the  tracking  of  multiyear  ice  features  to  infer  ice 
motion. 

No  evidence  of  altered  ice  morphology  in  the  Yermak  Plateau  region  was 
detected.  Either  the  reported  anomalous  currents  do  not  affect  ice  morphology  or 
the  limitations  inherent  to  passive  microwave  imagery  (poor  spatial  resolution) 
prevent  detection.  On  the  other  hand,  distinct  multiyear  ice  features  are  observed 
to  move  in  the  direction  of  mean  ice  flow  and  are  well  correlated  with  the  motion 
of  nearby  data  buoys.  Short  term  ice  motion  vectors  are  unreliable  because  of  the 
low  spatial  resolution  of  the  passive  microwave  data.  However,  long  term 
averages  are  indicative  of  ice  flow  trajectories  and  can  be  used,  for  example,  to 
increase  our  knowledge  of  the  interannual  variability  of  heat  and  salinity  budgets 
in  specific  regions,  including  Fram  Strait,  the  principal  outflow  region  of  the  arctic 
basin. 
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I.    INTRODUCTION 

All-weather,  continuous  observations  of  the  arctic  ice  cover  are  important  for 
commercial  and  naval  operations  and  for  scientific  studies  such  as  global  climate 
change.  However,  such  observations  are  generally  difficult  to  obtain.  In-situ 
measurements  from  ice  camps  and  ships  are  expensive  to  acquire  and  spatially 
limited.  Aircraft  measurements  can  expand  the  areal  coverage,  but  asset  and 
resource  limitations  prevent  continuous  large-scale  monitoring  of  the  ice  pack. 
Satellite  platforms  are  the  only  practical  means  of  providing  continuous  coverage. 
Of  the  passive  sensing  methods,  only  microwave  emissions  penetrate  the  clouds 
and  darkness  that  shield  much  of  the  arctic  ice  from  visual  and  infrared 
observations.  Passive  satellite  microwave  data  are  operationally  available  from 
Special  Sensor  Microwave/Imager  (SSM/I)  sensors  on  the  Defense 
Meteorological  Satellite  Program  (DMSP)  series  satellites  (Hollinger  et  al.,  1987). 
Historical  data  are  also  available  to  researchers  on  CD-ROM  distributed  by  the 
National  Snow  and  Ice  Data  Center  (NSIDC). 

The  objectives  of  this  study  are  to  explore  the  use  of  passive  microwave 
satellite  measurements  for  arctic  sea  ice  studies.    Two  applications  are  examined: 

•  The  detection  of  altered  ice  morphology  due  to  anomalous 
diurnal  tidal  currents  observed  over  the  Yermak  Plateau.  If 
successful,  this  technique  can  be  applied  in  other  polar  regions  to 
expand  our  knowledge  of  the  local  tidal  structure  and  its  effect  on 
sea  ice  cover. 

•  Determination  of  ice  flow  by  tracking  the  movement  of  floes 
exhibiting  different  morphologies  from  the  surrounding  ice  field. 


This  can  augment  the  Arctic  Buoy  Program  and  increase  our 

understanding  of  ice  dynamics. 
Chapter  II  presents  background  material  on  the  SSM/I,  data  sources  and  sea 
ice  microwave  radiation  characteristics.  Determination  of  sea  ice  concentrations 
from  passive  microwave  data  are  also  discussed.  Chapter  III  examines  the 
anomalous  diurnal  currents  observed  over  the  Yermak  Plateau.  Previous 
observations  and  models  prove  the  existence  and  explain  the  physics  behind 
these  currents.  The  effect  of  these  currents  on  the  ice  cover  is  examined  through 
passive  microwave  emissions.  Chapter  IV  shows  how  ice  features  can  be  tracked 
to  infer  general  ice  drift  patterns  in  the  arctic.  When  combined  with  other  data, 
this  technique  can  help  researchers  understand  ice  and  salinity  budgets  in  the 
polar  region.  Chapter  V  summarizes  the  results  and  closes  with  a  discussion  of 
further  research  that  could  be  attempted  using  these  data. 


II.  BACKGROUND 

A.  THE  SSM/I 

The  SSM/I  is  flown  on  the  Defense  Meteorological  Satellite  Program  (DMSP) 
block  5D-2  F8  satellite,  launched  on  June  19,  1987.  It  is  in  circular  sun- 
synchronous  near-polar  orbit  at  an  altitude  of  approximately  833  km,  an 
inclination  of  98.8°  and  a  period  of  102  minutes.  This  produces  14.1  full  orbit 
revolutions  per  day  with  a  0613  local  ascending  node  equatorial  crossing.  The 
SSM/I's  conical  scan  covers  a  1400  km  swath  which  provides  24  hour  near-global 
coverage.  The  orbit  inclination  precludes  coverage  of  a  280  km  radius  region 
centered  on  each  pole  extending  to  about  87.5°N.  Equatorial  regions  require  72 
hours  for  full  coverage,  however  overlapping  orbits  provide  full  daily  coverage 
poleward  of  45°N.  Hollinger  et  al.  (1990)  provide  a  good  discussion  of  SSM/I 
operating  characteristics. 

The  SSM/I  uses  four  frequencies  to  measure  seven  channels  of  linearly 
polarized  radiation  as  shown  in  Table  1.  Footprint  resolution  increases  with 
increasing  frequency.  Unfortunately,  the  high  resolution  85  GHz  vertically 
polarized  (85V)  channel  was  damaged  when  the  satellite  overheated  in  the  winter 
of  1987.  The  channel  was  useless  by  the  end  of  January  1989  (Hollinger  et  al., 
1990).  The  National  Snow  and  Ice  Data  Center  (NSIDC)  recently  reported  that 
the  85  GHz  horizontally  polarized  (85H)  channel  is  becoming  increasingly 
unreliable  and  is  basically  useless  (Taylor,  1990).  Poe  and  Conway  (1990) 
reported  the  current  SSM/I  geo-location  error  to  be  approximately  10-15  km. 


NSIDC  applies  a  correction  algorithm  that  reduces  this  error  to  8  km  (NSIDC, 

1990). 


TABLE  1:    SSM/I  FREQUENCIES 


Channel 

Frequency 
(GHz) 

Wavelength 
(cm) 

Polarization 

Resolution 
(km) 

1 

19.35 

1.55 

Vertical 

56 

2 

19.35 

1.55 

Horizontal 

56 

3 

22.235 

1.35 

Vertical 

45 

4 

37.0 

0.81 

Vertical 

33 

5 

37.0 

0.81 

Horizontal 

33 

6 

85.5 

0.35 

Vertical 

14 

7 

85.5 

0.35 

Horizontal 

14 

SSM/I  data  are  transmitted  to  Telemetry  and  Data  Relay  Stations  on  each  orbit. 
These  data  are  then  relayed  via  satellite  to  Navy  and  Air  Force  processing  centers 
which  produce  near  real  time  global  maps  of  sea  ice  location,  age,  type  and  other 
environmental  parameters.  Raw  data  are  also  continually  broadcast  to  Air  Force 
and  Navy  tactical  units  positioned  within  the  orbital  swath  of  the  satellite  that  are 
equipped  with  SMQ-1 1  satellite  receivers. 


Data  record  files  are  created  from  raw  radiometer  measurements.  Temperature 
data  records  (TDRs)  are  antenna  temperatures  which  are  calibrated  each  scan 
using  a  300°  K  black-body  radiator  and  the  3°  K  cosmic  radiation  as  reference 
temperatures.  An  antenna  pattern  correction  is  applied  to  the  TDRs  to  correct  for 
energy  received  in  the  sidelobes,  cross-polarization  coupling  and  feedhorn 
spillover  loss  (NSIDC,  1990).  These  main  beam  brightness  temperatures  are 
placed  in  sensor  data  record  (SDR)  files.  Environmental  parameters  are 
determined  by  applying  appropriate  algorithms  to  the  SDRs  and  are  stored  in 
environmental  data  records  (EDRs).  The  primary  archiver  of  SSM/I  data  in  orbital 
format  is  the  NOAA/National  Environmental  Satellite,  Data,  and  Information 
Service,  Satellite  Data  Services  Division  (NOAA/NESDIS/SDSD).  TDRs,  SDRs 
and  EDRs  are  also  forwarded  to  the  Naval  Research  Laboratory  for  archival. 

B.  SEA  ICE  RADIATION  CHARACTERISTICS 

Sea  ice  is  a  complex  material.  A  good  comprehensive  overview  of  both  its 
physical  and  microwave  radiation  characteristics  can  be  found  in  Parkinson  et  al. 
(1987).  The  following  is  a  general  discussion  of  these  properties. 

1.     Physical   Properties 

Sea  ice  covers  approximately  seven  percent  of  the  word's  oceans.  It  is 
inhomogeneous,  anisotropic  and  never  in  equilibrium.  Sea  ice  can  simultaneously 
contain  all  three  phases  of  matter:  solid  ice,  liquid  brine  pockets  and  air  cavities. 
The  fractions  of  the  constituent  phases  change  over  time.  New  ice  is  thinner  and 
saltier  than  first  year  ice,  which  transforms  after  a  melt  cycle  (summer  season)  into 
fresher,  thicker  and  less  dense  multiyear  ice  containing  numerous  air  pockets.  Ice 
can  be  covered  with  a  layer  of  dry  or  wet  snow  of  varying  thickness;  summer 


melting  forms  ponds  of  liquid  water  over  the  ice  sheet.    Figure  1  shows  the 
different  physical  characteristics  of  first  and  multiyear  ice. 

The  ice  composition  is  initially  determined  by  environmental  factors  at 
the  time  it  is  formed.  These  factors  include  air  temperature,  wind  speed,  ocean 
currents,  near  surface  salinity  and  weather  (Weeks,  1981).  The  ice  composition 
changes  continuously  due  to  weather,  ice  advection  and  fracturing  or  ridging 
effects.  Warmer  summer  temperatures  are  especially  important  as  melt  water 
percolates  through  the  ice  and  brine  pockets  are  flushed. 
2.     Radiative  Properties 

The  radiative  properties  of  sea  ice  are  determined  by  the  emissivity 
differences  between  ice  and  water.  Emissivity  is  the  ratio  of  emittance  (J  m~2  s_1) 
per  unit  wavelength  to  that  of  a  blackbody  at  the  same  temperature.  The 
Rayleigh-Jeans  approximation  (linear  relationship  between  emittance  and 
temperature)  applies  at  microwave  wavelengths.  The  emissivity  (e)  of  a  substance 

can  be  expressed  as: 


e=^  « 


where  Te  is  the  radiated  microwave  brightness  temperature,  proportional  to  the 
emittance,  and  Tsfc  is  the  physical  temperature  of  the  substance.  The  SSM/I 
radiometer  senses  Tb,  the  brightness  temperature  at  the  satellite  which  is  related 
to  Te  after  attenuation  and  contamination  by  atmospheric  effects  are  considered. 

Brightness  temperatures  can  be  used  to  calculate  ice  concentrations  and 
differentiate  between  open  water,  first  year  ice  and  multiyear  ice  based  on 
emissivity  variations  of  the  substances.    To  a  lesser  extent,  snow  and  water- 


covered  ice  can  also  be  observed.  Figure  2  summarizes  the  emissivity  differences 
between  sea  water,  first  year  ice  and  multiyear  ice  as  a  function  of  polarized 
frequency. 

The  ice  composition  controls  the  dielectric  properties  that  influence  the 
emissivity  of  the  ice.  Following  Comiso  (1983),  the  principal  difference  between 
sea  ice  and  freshwater  ice  are  brine  inclusions  in  the  former.  These  brine  pockets 
are  small,  typically  3-5  mm  by  0.25  mm  (Vant  et  al.,  1978).  The  brine  content  of 
sea  ice  depends  on  freezing  conditions  and  subsequent  temperature  changes  and 
varies  from  16  psu  for  first  year  ice  to  less  than  1  psu  for  multiyear  ice.  The  liquid 
brine  effects  the  dielectric  property  that  determines  emissivity  and  the  skin  depth, 
which  is  the  depth  from  which  most  of  the  observed  radiation  emanates.  In  saline 
ice  the  skin  depth  is  on  the  order  of  the  radiation  wavelength.  For  desalinated  ice 
this  depth  can  be  many  wavelengths  (Gloersen  et  al.,  1973). 

The  fundamental  difference  between  first  year  and  multiyear  ice  is  the 
presence  of  brine  pockets  in  the  former  and  their  replacement  in  the  freeboard 
layer  by  air  in  the  latter.  Because  of  the  small  skin  depth,  microwave  emissions 
from  the  more  saline  first  year  ice  emanate  from  the  snow/ice  interface;  emissions 
from  fresher  multiyear  ice  emanate  from  a  thicker  ice  layer  and  are  subject  to 
scattering  by  the  air  pockets.  This  scattering  reduces  the  net  radiation  from  the 
ice  surface  thus  lowering  the  brightness  temperature  (Comiso,  1983). 

Snow  cover  also  effects  microwave  radiation.  Snow  insulates  the  ice 
from  the  atmosphere,  leading  to  warmer  ice  from  oceanic  heat  conduction  and 
therefore  slightly  higher  brightness  temperatures.  This  may  be  compensated  by 
the  scattering  effect  of  snow  ice  crystals,  especially  at  higher  frequencies. 
Emissions  from  the  crystals  are  negligible  (Comiso  et  al.,  1982).    During  the 


summer  melt  season,  however,  each  crystal  is  surrounded  by  a  thin  layer  of  water. 
This  wet  snow  cover  becomes  the  primary  radiation  source,  significantly  lowering 
brightness  temperatures  due  to  water's  low  emissivity.  This  effect  makes  it  nearly 
impossible  to  discriminate  wet  snow-covered  multiyear  ice  from  first  year  ice  in 
the  summer. 

Ice  thickness  also  influences  brightness  temperatures.  At  microwave 
frequencies,  thin  ice  (radiation  wavelength  or  less)  is  contaminated  by  the  coldly 
radiating  water  below,  lowering  the  brightness  temperature.  This  may  be 
compensated  by  the  higher  salinity  and  warmer  physical  temperature  of  the  thin 
ice  versus  the  relatively  fresher  and  colder  thick  ice. 

Melt  ponds  are  formed  by  solar  heating  of  the  ice  cover  and  are 
characterized  by  thin  layers  of  water  scattered  among  solid  ice  regions.  The  size 
and  depth  of  the  ponds  depend  on  the  incident  solar  flux,  ice  salinity  and  ice 
roughness  (Goroch  and  Fett,  1990)  If  the  thin  water  layer  is  thicker  than  the  skin 
depth  of  the  radiation  (typically  a  few  mm  to  one  cm),  the  radiation  is  emitted 
primarily  by  the  liquid  and  has  a  much  colder  brightness  temperature  than  the  ice. 
For  thinner  layers  the  underlying  ice  will  contribute  some  radiation,  raising  the 
brightness  temperature.  The  variation  of  brightness  temperature  depends  on  the 
layer  depth  and  the  emission  frequency.  For  depths  on  the  order  of  a  few 
millimeters  the  19  and  22  GHz  channels  will  be  unaffected  (Goroch  and  Fett, 
1990).  As  the  frequency  increases,  contamination  by  the  water  layer  also 
increases.  The  total  effect  on  brightness  temperature  is  a  function  of  areal  extent 
within  the  radiometer  footprint  and  melt  pond  depth. 

Other  factors  affecting  microwave  emissivities  include  the  physical 
temperature  of  ice  and  undulations  caused  by  ridging  and  refreezing  of  melt 


ponds.   As  with  liquid  melt  ponds,  these  conditions  vary  spatially  and  temporally 
and  could  cause  observable  changes  in  ice  emissivity. 
3.     Atmospheric   Effects 

Important  atmospheric  constituents  which  affect  microwave  radiation 
are  water  vapor,  oxygen  and  liquid  water  content.  Seasonal  climatological 
constants  developed  from  physical  models  or  observations  are  generally  used  in 
algorithms  to  quantify  their  cumulative  effect  on  brightness  temperatures. 

Interaction  with  non-precipitating  clouds  in  the  arctic  region  is 
negligible  due  to  the  small  size  of  the  droplets  compared  with  the  microwave 
wavelengths  (Comiso,  1983).  Liquid  precipitation,  however,  can  contaminate 
brightness  temperatures.  Operational  algorithms  usually  flag  rain-contaminated 
radiances  for  rejection. 

Although  the  arctic  regions  have  low  humidities,  absorption  by  water 
vapor  can  effect  remotely  sensed  emissions.  For  the  SSM/I  only  the  22  and  85.5 
GHz  channels  are  significantly  effected  by  water  vapor  absorption  (Hollinger  et 
al,  1987).  Neither  channel  is  currently  used  in  ice  studies,  although  the  higher 
resolution  offered  by  operational  85.5  GHz  channels  in  future  DMSP  satellites 
will  be  tempting. 

The  atmospheric  oxygen  distribution  is  well  known  and  constant. 
Absorption  effects  are  included  in  seasonal  climatological  constants  used  in  ice 
concentration  algorithms. 

C.  SEA  ICE  CONCENTRATION  ALGORITHM 

1 .     Physical  Model 

A  simple  radiative  model  expressing  the  radiation  sensed  at  satellite 
height  (Tg)  is  presented  in  Svendsen  et  al.  (1983).    Surface  and  atmospheric 


effects  are  included,  plus  a  minor  free  space  contribution.  The  model  is  composed 
of  four  terms  (see  Figure  3):  radiation  from  the  surface  as  seen  through  the 
atmosphere  (A);  upwelled  radiation  from  the  atmosphere  (B);  downwelled 
atmospheric  radiation  reflected  at  the  surface  (C);  and  free  space  radiation 
reflected  by  the  surface  (D). 

TB  =  esfcTsfcO-^a)  +  STaxa  +  (l-esfc)5Taxa(l-xa)  +  (l-esfc)Tsp        (2) 
(A)  (B)  (C)  (D) 

Ta  is  the  optical  depth,  TSp  is  the  free  space  temperature  and  5Ta  is  the  weighted 
average  atmospheric  temperature  in  the  lower  troposphere.  Since  the  optical 
depth  is  small,  the  extinction  exponential  (e~Ta)  is  approximated  by  l-xa. 
Atmospheric  absorption  of  free  space  radiation  is  neglected. 

From  equation  (1)  the  emitted  brightness  temperature  (Te)  is  dependent 
only  on  the  physical  temperature  of  the  emitter  (Tsfc)  and  its  emissivity  (£sfc).  By 
substituting  Te/Tsfc  for  esfc  in  equation  (2),  Te  can  be  computed  if  xa  and  8Ta  are 
known,  and  if  Tsfc  is  inferred  from  5Ta. 

Three  components  contribute  to  Te:  sea  water,  first  year  ice  and 
multiyear  ice.  The  sum  of  the  fractional  components  of  each  (Cw,  Cf  and  Cm> 
respectively)  add  up  to  1 : 

l=Cm+Cf+Cw  (3) 

The  emitted  brightness  temperature  is  the  sum  of  the  individual  components, 
weighted  by  the  concentrations: 
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Te  =  £sfcTsfc  =  Cm^ice  +  Cf^ce  +  Cw£w(272)  (4) 

where  water  temperature  is  assumed  to  be  272°  K  and  £f,  em  and  ew  are  the 

emissivities  of  first  year  ice,  multiyear  ice  and  seawater,  respectively.    Tice  in 

equation  (4)  is  determined  by  assuming  a  linear  relationship  between  the 
atmospheric  surface  temperature  (Ta)  and  the  water  temperature: 

Tice  =  aTa  +  (l-a)(272)  (5) 

where  a  is  an  observationally  determined  constant  (=  0.4  from  Svendsen  et  al., 
1983). 

Fractional  concentrations  can  be  computed  from  equations  (1)  -  (5) 
using  two  radiometer  channels.  Sources  of  error  include  approximations  of 
optical  depth,  ice  temperature,  atmospheric  temperature  weighting  factor, 
atmospheric  surface  temperatures  and  assuming  constant  emissivities  for  first  and 
multiyear  ice.  Svendsen  et  al.  (1983)  compared  aircraft  radiometric  measurements 
with  Nimbus  7  radiometer  measurements  using  this  technique  with  in-situ 
temperature  observations  and  determined  accuracies  to  be  ±  3%  of  total  ice 
concentration  and  ±  10%  for  multiyear  ice  concentration  with  a  spatial  resolution 
of  90  km. 

2.     Multichannel   Algorithm 

a.    Multispectral   Analysis 

Comiso  (1983)  applied  a  multi-spectral  technique  that  exploits  the 
variation  of  emissivity  with  frequency  for  sea  water  and  various  ice  types  to 
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distinguish  between  various  first  and  multiyear  ice  and  open  water.  Figure  4 
shows  three  distinct  brightness  temperature  clusters  that,  when  geo-located, 
correspond  to  first  year  ice,  a  first/multiyear  ice  mixture  and  open  ocean 
conditions.  The  line  that  connects  the  first  year  cluster  with  the  open  ocean  is 
interpreted  as  varying  ice  concentrations.  Because  the  skin  depth  is  frequency 
dependent,  the  grouping  beneath  the  first  year  ice  is  theorized  to  be  thin  ice, 
where  cold  contamination  from  sea  water  is  observed  in  the  lower  (deeper 
penetrating)  frequency.  The  larger  spread  of  emissivities  and  brightness 
temperatures  in  the  37V  channel  reflects  the  greater  sensitivity  of  37  GHz 
radiation  to  near-surface  effects  due  to  its  smaller  skin  depth. 
b.   Polarization/Gradient   Ratios 

Cavelieri  et  al.  (1984)  used  polarization  ratios  (PR)  for  a  single 
frequency  and  spectral  gradient  ratios  (GR)  between  two  frequencies  to  estimate 
ice  concentrations  and  fractions.  PR  and  GR  are  defined  as  follows: 


PR"TBv(fl)  +  TBh(fl)  () 


TBv(fi)  +  TBv(f2)  W 


where  TBv  and  TBh  are  brightness  temperatures  for  vertical  and  horizontal 
polarizations  respectively  and  fj  and  f2  are  separate  frequencies.  Table  2  shows 

PRs  and  GRs  for  sea  ice  and  water  emissivities  at  the  two  SSM/I  frequencies  (37.0 
and  19.35  GHz)  used  to  calculate  ice  concentrations  and  fractions.  Note  that  PR 
is  much  greater  over  the  open  ocean  than  over  sea  ice.    Thus  PR  is  a  good 


12 


water/ice  discriminator  and  can  be  used  to  estimate  ice  concentrations.  GR 
depends  on  ice  type  and  is  affected  by  season.  Ice  types  are  determined  by  the 
sign  of  GR:  positive  over  the  open  ocean,  negative  over  multiyear  ice  and  near 
zero  over  first  year  ice  (Parkinson  and  Cavalieri,  1989). 


TABLE  2:  EMISSIVITIES  AT  19.35  AND  37  GHZ 


Ice  Type 

First  Year 
Thin  Ice 

First  Year 
Ice 

Multiyear 
Ice 

Calm 
Sea  Water 

Polarization 

V 

H 

V 

H 

V 

H 

V 

H 

Frequency 
(GHz) 

■ 

19.35 

0.96 

0.78 

0.97 

0.84 

0.86 

0.73 

0.60 

0.33 

37.0 

0.96 

0.87 

0.97 

0.95 

0.69 

0.64 

0.70 

0.40 

Polarization 
Ratio 

(37V  -  37H) 

0.09 

0.02 

0.05 

0.30 

Gradient 
Ratio 

(37V- 19V) 

0.00 

0.00 

-0.17 

0.10 

Cavalieri  et  al.  (1984)  describe  how  these  ratios  are  used  to 
calculate  total  and  multiyear  ice  concentrations  and  the  physical  temperature  of 
the  sea  surface.  They  conclude  that  the  precision  of  the  calculated  ice 
concentrations  ranges  from  5  -  9  %  and  the  precision  of  the  calculated  multiyear 
ice  fractions  vary  from  13  -  25  %.  A  caveat  on  the  ice  fraction  calculation  is  that 
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first  and  multiyear  ice  can  not  be  discriminated  during  the  summer  season 
because  the  general  surface  melt  conditions  mask  the  distinguishing  properties  of 
the  two  ice  types.  A  non-iterative  version  of  this  method  has  been  adopted  by 
the  NASA  Sea  Ice  Algorithm  Working  Group  and  is  used  by  the  NSIDC  (NSIDC, 

1990). 

D.  DATA  SOURCES 
1.     NSIDC  CD-ROMs 

NSIDC  distributes  SSM/I  brightness  temperatures  on  CD-ROM  disks  for 
both  arctic  and  antarctic  regions.  A  continuous  record  is  currently  available  from 
9  July  1987  -  31  December  1989.  Earlier  SMMR  data  are  available  in  the  same 
format  from  1978-1987.  Subsequent  SSM/I  data  are  in  production  for  CD-ROM 
distribution.  The  data  include  all  SSM/I  channels  with  85  GHz  data  gridded  at 
12.5  km  resolution  and  19,  22  and  37  GHz  data  gridded  at  25  km  resolution. 
Added  features  are  land  masks  and  coast  outlines  that  provide  geographic 
reference  points  on  the  imagery.  The  data  are  obtained  from 
NOAA/NESDIS/SDSD.  Fleet  Numerical  Oceanography  Center  produces  the 
original  TDR  files  from  the  raw  satellite  data  as  part  of  their  operational  data 
processing.  The  TDRs  are  transferred  to  NOAA/NESDIS/SDSD  and  then  to 
NSIDC.  At  NSIDC  the  TDRs  are  examined  for  bad  data  values  and  a  geo- 
location  correction  is  applied.  The  antenna  temperatures  are  then  converted  to 
brightness  temperatures. 

Swath  data  from  each  channel  are  mapped  onto  the  appropriate  grid  by 
converting  the  geodetic  latitudes  and  longitudes  at  the  center  of  each 
observation  footprint  into  map  grid  coordinates.  Brightness  temperatures  are 
summed  in  each  grid  cell  over  a  24  hour  period  and  then  divided  by  the  number 
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of  observations.    The  resulting  data  fields  are  daily  averaged  mosaics  of  each 
polar  region. 

The  data  are  mapped  onto  polar  stereographic  projections  with  the  grid 
tangent  to  the  earth  at  70°  latitude  rather  than  at  the  poles.  This  increases  the 
distortion  at  the  poles  by  3%  but  is  offset  by  nearly  true  projections  in  the 
marginal  ice  zones.  Figure  5  shows  the  area  of  coverage  for  the  Northern 
Hemisphere. 

2.     Ice  Buoy  Tracks 

The  Arctic  Buoy  Program  has  been  operating  since  1979.  Objectives  of 
the  program  are  to  provide  surface  atmospheric  pressure  measurements  over  the 
arctic  basin  and  to  define  the  large  scale  field  of  sea  ice  motion  (Thorndike  and 
Colony,  1980).  Buoys  are  parachuted  from  C-130  aircraft  at  designated 
locations.  Each  buoy  transmits  pressure  and  air  temperature  information  to  a 
satellite  at  one  minute  intervals.  This  information  is  only  received  when  the 
satellite  is  within  the  radio  field  of  view.  The  data  are  retransmitted  to  a  ground 
station,  decoded  and  filtered,  and  forwarded  to  the  Polar  Science  Center  at  the 
University  of  Washington.  Thorndike  and  Colony  (1980)  examined  the 
locational  record  of  one  buoy  and  reported  the  positional  accuracy,  as  measured 
through  the  Doppler  shift  of  the  buoy  signal  received  by  the  satellite,  to  be  within 
300  m  for  two  thirds  of  the  reported  positions.  Very  few  fixes  were  reported  with 
errors  greater  than  10  km  and  these  can  be  easily  removed  from  the  record  by 
careful  editing. 
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Figure  1  Schematic  illustration  of  differences  between  (a)  first  year  ice  and 

(b)  multiyear  ice  (after  Comiso,  1983). 
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Figure  2  Emissivities  of  calm  seawater,  multiyear  ice  and  First  year  ice  as  a 
function  of  frequency  for  horizontal  and  vertical  polarizations 
(modified  from  Svendsen  et  al.,  1983). 
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Figure  3         Illustration  of  the  four  terms  in  the  radiation  model. 
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Figure  4         Multispectral  analysis  of  brightness  temperatures  (from  Comiso, 
1983). 
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Figure  5         SSM/I  north  polar  grid  (from  NSIDC,  1990). 
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HI.      CASE  STUDY:  ANOMALOUS  DIURNAL  CURRENTS  ON 

THE  YERMAK  PLATEAU 

A.   BACKGROUND 

1.  The  Yermak  Plateau 

The  Yermak  Plateau  (Figure  6)  is  a  hooked,  finger-like  extension  of  the 
Barents  Sea  continental  slope,  located  off  the  northwestern  coast  of  Svalbard. 
The  plateau  extends  400  km  into  the  Arctic  Ocean  and  covers  roughly  20,000 
km^.  Plateau  depths  range  from  1000  to  600  m.  The  shallowest  depth  is  425  m 
on  a  ridge  striking  north-northwest  near  the  plateau's  central  eastern  edge  (Perry, 
1986).  The  plateau  slopes  down  rapidly  to  3500  m  to  the  north  and  northwest. 
The  eastern  slopes  gradually  level  out  to  a  semi-enclosed  2000  m  deep  basin. 

2.  Regional  Tides 

Tides  in  the  Arctic  Ocean  are  predominantly  semi-diurnal.  Schwiderski 
(1986)  developed  a  model  for  global  tidal  prediction  based  on  a  1°  by  1°  grid 
with  accuracies  better  than  10  cm.  Figures  7  and  8  show  the  structure  of  the 
leading  semi-diurnal  (M2)  and  leading  diurnal  (Ki)  partial  tides,  respectively,  for 
the  Nordic  Seas.  Dashed  lines  depict  the  co-range  lines  measured  in  centimeters. 
The  M2  amplitude  varies  from  zero  at  the  amphidromic  points  to  80  cm  and 
ranges  from  20-30  cm  over  the  Yermak  Plateau.  The  amplitude  of  the  Ki  tide 
varies  only  from  zero  to  10  cm,  ranging  from  2-5  cm  over  the  Plateau. 

Gjevik  and  Straume  (1989)  developed  a  regional  tidal  model  to  simulate 
the  M2  and  Ki  tides  for  the  Arctic  Ocean  and  Nordic  Seas  with  a  50  by  50  km 
grid  size.   The  M2  and  Kl  amplitudes  are  similar  to  Schwiderski's  global  model. 
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The  regional  model  also  shows  clockwise  Kl  current  rotation  along  the  west 
Norwegian  shelf  up  to  Svalbard.     These  authors  believe  this  indicates  the 
presence  of  a  coastal  trapped  wave  modal  structure. 
3.     Observations 

Anomalous  diurnal  currents  over  the  Yermak  Plateau  were  first  reported 
by  Hunkins  (1986)  from  FRAM  III  and  IV  current  meter  data.  The  meters  drifted 
from  deep  (4000  m)  water  to  shallow  water  over  the  plateau's  northwestern  flank. 
Hunkins  resolved  absolute  currents  into  orthogonal  components  in  the  along 
slope  (60°)  and  cross  slope  (150°)  directions  and  found  that  strong  diurnal 
currents  exist  along  the  plateau's  slope.  Figure  9  shows  current  velocities  at  104 
m  along  with  the  bathymetric  profile  for  the  FRAM  IV  drift.  During  the  deep 
water  portion  of  the  drift,  the  current  record  is  noisy  with  some  evidence  of  semi- 
diurnal tidal  motion.  Over  the  slope  diurnal  motion  dominates  the  signal.  Once 
on  top  of  the  plateau,  the  semi-diurnal  signal  returns.  Kinetic  energy  spectra 
confirm  the  dominance  of  diurnal  over  semi-diurnal  amplitudes  (Figure  10). 

Tidal  motion  polarization  differs  for  the  two  tidal  periods.  Figure  1 1 
shows  a  broad  diurnal  peak  in  the  kinetic  energy  spectra  near  150°  and  a 
minimum  near  60°.  Conversely,  the  semi-diurnal  motions  are  nearly  parallel  to  the 
slope  with  a  maximum  at  60°  and  a  minimum  at  150°.  The  semi-diurnal  alignment 
agrees  with  the  modeled  M2  tides  for  the  Arctic  Ocean  (Kowlik,  1979  as  cited  by 

Hunkins,  1986).    Hunkins  believed  the  cross-shelf  diurnal  currents  indicated 
strong  topographic  vorticity  waves. 

Strong  diurnal  currents  were  also  reported  from  MIZEX  84  data. 
Morison  et  al.  (1987)  noted  semi-diurnal  currents  dominated  over  most  of  the 
current  meter  drift  record  with  amplitudes  of  10  cm/s.  Diurnal  motions,  however, 
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predominated  when  the  current  meter  drifted  over  the  southwestern  slope  of  the 
Yermak  Plateau.  The  diurnal  component  of  motion  ranged  up  to  12  cm/s. 

The  observations  of  Hunkins  (1986)  and  Morison  et  al.  (1987)  were 
confirmed  by  Hoffman's  (1990)  analysis  of  ARCTEMIZ  86  drifter  buoy  tracks. 
He  showed  a  strong  correlation  between  the  enhanced  diurnal  currents  and 
bathymetry.  While  in  deep  water  no  significant  peak  in  the  energy  spectrum  at 
the  semi-diurnal  or  diurnal  period  was  found.  However,  once  the  buoys  crossed 
the  3000  m  contour  into  shallower  water  the  energy  spectra  exhibited  a  diurnal 
peak  which  exceeded  the  semi-diurnal  peak.  Figure  12  is  a  buoy  track  that 
clearly  shows  diurnal  clockwise  looping  over  the  plateau's  northern  slope. 
Hoffman  also  reported  a  periodic  increase  in  diurnal  clockwise  motion  that  he 
associated  with  the  fortnightly  (spring-neap)  beat  in  surface  height  inequality 
caused  by  the  superposition  of  the  two  principal  semi-diurnal  tidal  components. 
4.     Analysis 

Hunkins  (1986)  showed  that  topographically  trapped  vorticity  waves 
can  exist  at  the  diurnal  frequency  over  the  Yermak  Plateau.  He  ignored  the 
details  of  forcing  and  used  a  simple  free  wave  analytic  model.  Chapman  (1989) 
explored  the  actual  forcing  of  the  waves  by  the  diurnal  tide.  He  showed  that  a 
rectilinear  oscillating  barotropic  current  is  significantly  amplified  by  Hunkin's 
idealized  representation  of  Yermak  Plateau  topography. 

Recent  work  by  Prazuck  (1991)  indicated  high  kinetic  energy  cells  are 
possible  on  the  Plateau.  His  linear  barotropic  numerical  model  created  such  cells 
in  close  agreement  with  observed  diurnal  currents.  He  showed  through  an 
analytic  model  that  the  highest  energy  cell  could  be  due  to  the  backscattered 
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energy  of  a  coastally  trapped  wave  that  is  unable  to  propagate  beyond  the  tip  of 
the  plateau. 

Prazuck  concluded  that  the  enhanced  rotating  diurnal  currents  may 
influence  ice  floe  motion.  The  high  energy  cells  would  increase  the  grinding 
together  of  the  floes  and  alter  the  break-up  and  freeze-up  patterns  in  the  region. 
Several  phenomena  might  be  related  to  this  interaction: 

•  Strong  turbulent  fluxes  associated  with  the  diurnal  currents 
such  as  those  found  during  CEAREX  89  by  Padman  and  Dillon 
(1991). 

•  Enhanced  ambient  noise  in  the  region  of  grinding.  Cousins 
(1991)  correlated  strong  variations  in  ambient  noise  levels  with 
12  hour  tidal/inertial  oscillations  east  of  the  Yermak  Plateau. 

•  Altered  ice  morphology  due  to  the  convergence/divergence 
of  icefloes. 

5.     Effects  on  Sea  Ice  Characteristics 

The  interaction  of  the  enhanced  diurnal  currents  on  the  sea  ice  should 
be  a  seasonal  phenomena.  The  high  ice  concentrations  of  the  winter  pack  ice 
should  significantly  dampen  floe  motion,  reducing  the  possibility  of  altering  ice 
conditions.  Lower  ice  concentrations  in  summer  on  the  other  hand  should  permit 
larger  amplitude  motions  and  increase  the  chances  of  remotely  observable 
changes  of  sea  ice  characteristics. 

Ice  floe  motion  effects  both  ice  concentration  and  morphology.  The 
enhanced  currents  can  induce  local  divergence  in  the  ice  field  and  raise  or  lower 
ice  concentrations  on  a  scale  observable  from  space.  Divergence  and  grinding 
together  of  floes  can  increase  fracturing.    Subsequent  refreezing  would  increase 
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the  fraction  of  first  year  ice  in  a  predominantly  multiyear  ice  field  without 
significantly  affecting  the  total  ice  concentration. 

B.  METHODOLOGY 

Data  from  the  NSIDC  SSM/I  CD-ROM  data  set  were  used  to  test  the 
hypothesis  that  changes  in  ice  concentration/morphology  could  be  detected  by 
satellite  remote  sensing.  SSM/I  data  were  used  because,  although  visible  and  IR 
measurements  provide  higher  resolution  data,  cloud  contamination  impede  most 
attempts  to  examine  a  time  series  extending  beyond  a  few  days. 

Sea  ice  concentrations  were  computed  using  the  algorithm  adopted  by  the 
NASA  Sea  Ice  Algorithm  Working  Group  (NSIDC,  1990).  Gradient  ratios  were 
also  calculated  by  modifying  the  ice  concentration  program  to  return  gradient 
ratios  instead  of  ice  concentrations.  The  25  km  grid  spacing  used  by  the  NSIDC 
data  set  was  retained.  Appendix  A  contains  the  software  code  listing. 

Figure  13  shows  how  the  Yermak  Plateau  region  was  divided  into  test  and 
control  areas.  A  strip  of  5  pixels  (125  by  25  km),  extending  over  the 
northwestern  flank  of  the  plateau,  was  selected  as  the  test  area.  Control  areas  of 
the  same  size  were  located  50  km  to  the  northwest  and  125  km  to  the  southeast. 
These  strips  will  be  referred  to  as  'Slope',  'North'  and  'South',  respectively.  A 
perpendicular  strip  75  km  upstream  of  the  general  ice  flow,  northeast  of  the 
plateau  tip,  was  used  to  establish  the  gradient  ratio  of  ice  flowing  into  the  test 
area. 

A  168  day  time  series  of  daily  ice  concentrations,  extending  from  1  May  to  15 
October  1988,  was  studied.  One  four  day  and  several  one  to  two  day  data  gaps 
in  the  record  were  filled  by  linear  interpolation  of  surrounding  values.  As  seen  in 
Figure  14,  ice  concentrations  varied  from  75-100%  during  the  period  for  the 
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northern  and  test  strips,  and  dipped  to  as  low  as  35%  for  the  southern  strip.  At 
other  times  of  the  year  concentrations  remained  near  100%.  For  the  gradient  ratio 
analysis  the  time  series  extended  from  13  January  to  30  April  1988.  This  period 
avoided  brightness  temperatures  contaminated  by  melt  ponds. 

C.   RESULTS 

a.     Ice  Concentration 

The  ice  concentration  algorithm  was  first  checked  for  spatial 
consistency.  Figure  15  shows  concentration  values  for  each  pixel  within  the 
Slope  strip  during  the  period  of  interest.  Table  3  contains  a  correlation  matrix 
showing  that  the  values  are  spatially  consistent  with  0.74  as  the  lowest 
correlation  between  the  two  endpoint  pixels. 

Differencing  between  the  North  -  Slope  and  South  -  Slope 
concentration  values  (Figure  16)  showed  no  significant  deviations  between  ice 
concentrations  over  the  slope  and  control  areas.  The  North  -  Slope  mean 
difference  is  1.04%,  well  inside  the  10%  accuracy  of  the  concentration  algorithm. 
The  Slope  -  South  mean  difference  is  larger  at  8.56%,  but  as  seen  from  the 
concentration  record  this  is  due  to  the  generally  lower  concentrations  in  the 
South  strip  later  in  the  summer. 

The  record  was  next  examined  for  evidence  of  the  fortnightly  beat  of 
current  intensity  reported  by  Hoffman  (1986).  If  the  currents  influence  ice 
concentrations,  the  largest  signal  should  be  present  when  the  currents  are 
strongest.  An  autocorrelation  of  Slope  concentration  values  showed  no 
significant  correlations  over  a  fortnightly  period  (Figure  17).  Spectral  densities 
were  also  computed  for  the  mean  Slope,  North  and  South  concentrations,  as  well 
as  the  North  -  Slope  and  South  -  Slope  differences  (Figures  18-19).   In  none  of 
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these  cases  was  a  spectral  peak  located  at  or  near  the  13.66  day  lunar  fortnightly 
cycle  (0.07  cycles  per  day). 


TABLE  3:  CORRELATION  MATRIX  FOR  YERMAK 

SLOPE  PIXELS 


Position 

Pixel  1 

Pixel  2 

Pixel  3 

Pixel  4 

Pixel  5 

Pixel  1 

1 

Pixel  2 

.955 

1 

Pixel  3 

.881 

.954 

1 

Pixel  4 

.803 

.890 

.963 

1 

Pixel  5 

.743 

.823 

.915 

.966 

1 

b.     Gradient  Ratio 

The  gradient  ratio  is  indicative  of  ice  age.  As  shown  in  Table  2,  the 
gradient  ratio  should  be  near  zero  for  first  year  ice,  less  than  zero  for  multiyear  ice 
and  greater  than  zero  for  sea  water.  The  Siberian  Branch  of  the  Polar  Drift  Stream 
advects  mostly  multiyear  ice  through  this  region  (Polarstern  Scientific  Party, 
1988).  If  interactions  between  the  diurnal  currents  and  the  ice  are  sufficiently 
strong  to  increase  ice  fracturing  over  the  Yermak  Plateau  slope,  the  refreezing  of 
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the  resultant  small  leads  should  increase  the  total  fraction  of  first  year  ice.  An 
increased  fraction  of  more  saline  ice  will  drive  the  gradient  ratio  more  positive 
compared  with  the  control  areas  to  north,  south  and  upstream  of  the  flow. 

Gradient  ratios  for  the  North,  Slope  and  South  regions  are  shown  in 
Figure  20.  Mean  values  are  -0.033,  -0.031  and  -0.034,  respectively.  All  gradient 
ratios  are  negative  as  expected  from  the  predominantly  multiyear  ice.  Although 
the  gradient  ratios  are  at  least  26%  higher  over  the  slope  than  over  the  control 
areas  during  days  10  -  42,  the  difference  between  upstream  and  slope  values  is 
negligible  for  the  same  period  (Figure  21).  This  means  that  no  change  measurable 
by  the  SSM/I  occurred  as  the  ice  flowed  over  the  slope  region. 

D.   CONCLUSIONS 

Anomalous  diurnal  currents  over  the  Yermak  Plateau  slopes  have  been 
reported  and  analyzed  by  Hunkins  (1986),  Morison  et  al.  (1987),  Hoffman  (1990) 
and  Prazuck  (1991).  These  currents  have  the  potential  for  affecting  the  sea  ice 
cover  in  this  region  in  several  ways.  Divergence/convergence  patterns  could 
alter  local  ice  concentrations.  Stresses  induced  in  the  ice  could  lead  to  increased 
grinding  and  fracturing  of  floes,  changing  the  morphology  of  the  ice  cover.  If 
such  changes  could  be  noted  through  remote  sensing,  existing  data  bases  could 
be  explored  for  similar  phenomena  in  other  ice-covered  regions. 

Ice  concentration/morphology  changes  are  not  observed  by  the  SSM/I.  Two 
explanations  could  account  for  this  lack  of  confirmation: 

1)  The  currents  have  little  or  no  effect  on  the  ice  cover. 

2)  SSM/I  limitations  obscure  features  in  the  ice  caused  by  the  currents. 
Factors  supporting  the  second  hypothesis  include  poor  SSM/I  spatial  resolution 
and  contamination  by  melt  ponds  or  snow  cover.   Passive  microwave  emissions 
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observed  by  the  SSM/I  are  too  imprecise  to  verify  the  anomalous  current  effects 
on  the  sea  ice  cover. 
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Figure  6         Bathymetric  chart  of  the  Yermak  Plateau  (from  Hoffman,  1990). 
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Figure  7  Leading  semi-diurnal  (M2)  ocean  tide  (from  Schwiderski,  1986). 
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Figure  8         Leading  diurnal  (Ki)  ocean  tide  (from  Schwiderski,  1986). 
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Figure  9         Absolute  current  velocities  along  the  FRAM  IV  track  at  104  m  (from 
Hunkins,  1986). 
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Figure  10       Kinetic  energy  spectra  for  the  cross  slope  velocity  component  from 
FRAM  IV  at  303  m  (from  Hunkins,  1986). 
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Figure  11         Peak  amplitudes  of  kinetic  energy  for  the  Ki  and  M2  frequencies  as  a 
function  of  current  axis  orientation  for  FRAM  IV  at  303  m  (from 
Hunkins,  1986). 
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Figure  12        Diurnal  looping  of  ARCTEMIZ  86  drifter  buoys  (from  Hoffman, 
1990). 
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Figure  13        Yermak  Plateau  test  and  control  areas.     Hatched  boxes  indicate 
individual  pixels.  Depths  are  in  meters. 
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Figure  14       Total  ice  concentrations,  1  May  -  15  October  1988:  (a)  North 
control,  (b)  Slope  and  (c)  South  control. 
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Figure  15        Yermak  slope  ice  concentrations,  1  May  -  15  October  1988. 


39 


15. 


£   101 
5-1 


j  (a)  North  -  Slope  Difference 


o 

c 
o 
U 

o 


0-. 


■5-. 


-10 


B 


i— i—    i   i   i 


120   140   160   180   200   220   240   260   280   300 


IP- 


t 


50- 


45- 
40- 
35- 


304 


# 


o 

c 

o 
U 


25- 


204 


15 


(b)  Slope  -  South  Difference 
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Figure  16       Percent  concentration  differences  (1  May  -  15  October  1988):    (a) 
North  control  minus  Slope,  (b)  South  control  minus  Slope. 
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Figure  17        Autocorrelation  of  slope  region  ice  concentrations. 
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Figure  18       Ice  concentration  spectra  (1  May  -  15  October  1988):    (a)  North 
Control,  (b)  Slope,  and  (c)  South  control. 
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Figure  19        Ice  concentration  spectra  (1  May  -  15  October  1988):    (a)  North 
Slope  difference  and  (b)  Slope  -  South  difference. 
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Figure  20       Gradient  ratios  for  Slope  and  Control  areas,  13  January  -  30  April 
1988. 
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Figure  21        Slope  -  Upstream  gradient  ratio  difference,  1  February  -  30  April 
1988. 
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IV.     CASE  STUDY:     ICE  MOTION  ANALYSIS 

A.  BACKGROUND 

1 .     General  Ice  Flow 

Sea  ice  motion  affects  a  variety  of  important  arctic  processes.  Carsey 
and  Holt  (1987)  note  that  ice  movement  affects  latent  heat  advection,  oceanic 
surface  stress,  open  water  production  through  divergence  and  shear,  is  a  passive 
tracer  of  ocean  currents,  and  affects  icebreaker  navigation  and  ice  structure. 
Within  the  central  ice  pack,  relative  ice  motion  creates  both  open  water,  through 
which  significant  heat  is  lost,  and  ridges  that  are  operational  hazards.  At  the  ice 
margins  ice  growth  associated  with  divergence  affects  brine  production  which,  in 
turn,  influences  intermediate  and  deep  water  formation. 

Ice  motion  is  primarily  due  to  local  surface  wind  stress,  ocean  currents, 
Coriolis  and  inertial  forces,  and  effects  of  horizontal  stresses  transferred  through 
the  ice  itself  (Hibler,  1986).  Colony  and  Thorndike  (1984)  examined  a  90  year 
record  of  ship,  research  station  and  buoy  tracks  to  create  a  quantitative  estimate 
of  the  mean  ice  motion  in  the  arctic  basin  (Figure  22).  Two  features  dominate  the 
circulation  pattern — the  anticyclonic  Beaufort  Gyre  in  the  Canadian  Basin  and 
the  Transpolar  Drift  Stream  (TDS)  extending  across  the  Eurasian  Basin  to  Fram 
Strait. 

The  TDS  is  the  largest  and  most  concentrated  meridional  ice  flow  in  the 
world.  It  is  the  primary  outflow  region  for  arctic  sea  ice.  Variations  in  flow 
magnitude  are  theorized  to  effect  the  global  climate  (Vinje  and  Finnekasa,  1986). 
As  seen  in  Figure  23,  the  TDS  is  composed  of  two  branches.   The  Polar  Branch 
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advects  old  multiyear  ice  across  the  pole  and  out  through  Fram  Strait.  The 
Siberian  Branch  contains  thick,  deformed  multiyear  ice  floes  that  probably 
originated  on  the  Siberian  shelves.  Closer  to  the  Yermak  Plateau,  generally  thin, 
first  year  ice  formed  locally  or  in  the  Barents  and  Kara  Seas  is  found  (Polarstern 
Scientific  Party,  1988).  The  merging  of  these  different  ice  types  makes  Fram  Strait 
a  rich  area  for  SSM/I  feature  tracking. 
2.      Feature  Tracking 

Feature  tracking  estimates  motion  of  a  fluid  by  tracking  features  that  are 
advected  in  the  fluid.  Zhang  (1991)  identifies  two  main  steps  necessary  for 
successful  feature  tracking: 

1)  Tie  point  selection  -  Identification  in  the  first  image  of  suitable 
features  to  be  tracked;  and 

2)  Tie  point  registration  -  Matching  the  same  tie  points  in  a  subsequent 
image. 

The  underlying  assumption  is  that  the  tracked  features  are  advected  by 
the  fluid.  Motion  vectors  are  representative  of  the  surrounding  medium  only 
when  this  assumption  is  valid.  Possible  disruptive  influences  include  waves  in  the 
fluid,  when  fluctuations  in  fluid  properties  are  tracked  rather  than  the  Lagrangian 
flow  of  the  fluid,  atmospheric  contamination  which  can  temporarily  obscure 
features,  and  Eulerian  effects  such  as  local  heating  or  cooling  that  provide  feature 
sources  and  sinks. 

In  the  past  sea  ice  motion  has  been  inferred  from  the  drift  records  of  ice- 
locked  ships,  manned  research  stations  and  buoy  tracks.  The  number  and  spatial 
density  of  these  in-situ  observations  is  limited.    Satellites  offer  the  advantage  of 
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repeated  sampling  over  large  areas — a  potentially  rich  source  of  information  on 
the  spatial  and  temporal  structure  of  the  sea  ice  motion  field. 

Sea  ice  movement  can  be  computed  from  displacements  of  features  in 
sequential  images.  Subjective  techniques,  which  require  manual  tracking  of 
specific  floes,  are  labor  intensive  and  subject  to  repeatability  problems  due  to 
differences  in  operator  training.  Automated  techniques,  such  as  the  maximum 
cross  correlation  method  described  by  Ninnis  et  al.  (1986),  are  much  more  efficient 
but  must  be  heavily  filtered  to  remove  spurious  movement  vectors. 

Three  electromagnetic  regimes  are  currently  available  for  passive  space- 
based  observations  of  ice  motion.  In  order  of  decreasing  resolution,  they  are 
visible  (0.3  -  0.7  jim),  infrared  (3-12  (im)  and  microwave  (0.3  -  1.5  cm).  The  major 
limitations  of  high  resolution  visible  and  IR  imagery  are  their  inability  to  penetrate 
cloud  cover  and  the  light  dependence  of  visual  wavelengths.  For  example, 
Emery  et  al.  (1991)  investigated  Fram  Strait  ice  motion  with  IR  imagery  but  were 
limited  to  six  images  over  a  10  day  period  in  April  1986  that  contained 
exceptionally  cloud-free  conditions.  On  the  other  hand,  low  resolution  passive 
microwave  imagery  can  penetrate  clouds  and  are  independent  of  light  conditions. 

Active  microwave  systems  offer  the  advantages  of  cloud  penetration 
and  light  independence  with  much  higher  resolutions  than  passive  methods.  For 
example,  Synthetic  Aperture  Radar  (SAR)  data  from  SEASAT  had  a  resolution  of 
25  m  (Carsey  and  Zwally,  1986).  SAR  information  from  the  ERS-1  satellite  will 
soon  be  available,  however  its  application  will  be  limited  due  to  the  limited 
sampling  in  time  and  space,  especially  in  the  subpolar  regions  (Zhang,  1991). 

Many  researchers  believe  satellite  passive  microwave  sensing  is 
unsuitable  for  tracking  ice  motion  (Zhang,  1991;  Cavalieri,  1991   personal 
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communication)  due  to  its  inherent  coarse  spatial  resolution.  The  primary  use  of 
this  satellite  data  for  arctic  sea  ice  studies  has  been  to  examine  the  temporal  and 
spatial  variability  of  the  ice  cover.  Svendsen  et  al.  (1983)  reported  the  motion  of 
a  patch  of  high  multiyear  ice  concentration  that  correlated  well  with  buoy 
movements,  but  the  potential  for  feature  tracking  has  not  been  pursued  in  the 
published  literature  (Cavalieri,  1991  personal  communication). 

One  can  understand  why  passive  microwave  feature  tracking  has  been 
overlooked.  As  mentioned  earlier,  brightness  temperatures  are  sensitive  to 
variations  in  surface  conditions  such  as  wet  snow  or  melt  ponds.  The  spatial 
resolution  is  too  coarse  to  track  individual  floes.  Judicious  application  of  SSM/I 
data,  however,  can  overcome  these  limitations.  The  basic  requirements  are  large 
scale  spatial  gradients  in  ice  type  (as  determined  from  gradient  ratios)  to  overcome 
the  coarse  data  resolution,  and  small  temporal  variability  of  surface  conditions  to 
limit  melt  pond  or  wet  snow  contamination.  The  first  requirement  is  achieved  by 
examining  areas  where  old  and  new  ice  merge.  The  second  requirement  is  met  by 
rejecting  data  during  the  summer  melt  season. 

These  two  restrictions  limit  the  wide-scale  use  of  inferring  ice  motion 
through  SSM/I  feature  tracks.  Regions  where  this  technique  shows  the  most 
promise,  however,  are  areas  of  high  interest  such  as  Fram  Strait,  through  which 
occurs  approximately  90%  of  the  heat  exchange  and  75%  of  the  mass  exchange 
between  the  Arctic  and  World  Oceans  (Aagaard  and  Greisman,  1975). 

B.  METHODOLOGY 
1.     Satellite  Images 

Multiyear  ice  fractions  were  calculated  using  the  NASA  gradient  ratio 
technique  described  earlier.  The  same  software  that  computed  ice  concentrations 
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and  gradient  ratios  for  the  Yermak  Plateau  case  study  was  used  to  create 
multiyear  ice  fraction  images.  Appendix  A  contains  the  FORTRAN  code  listing. 
The  software  merges  the  computed  multiyear  ice  fractions  with  a  land  mask  to 
produce  a  304  by  448  pixel  polar  stereographic  representation  of  the  arctic 
region,  extending  from  the  pole  south  to  70°  N.  The  pixel  resolution  is  25  km. 

Images  were  produced  from  the  NSIDC  Winter  88  and  Winter  89  data 
sets.  Images  were  created  at  3  day  intervals.  If  one  assumes  a  typical  10  cm/s  ice 
speed,  a  feature  advected  by  the  flow  will  move  nearly  26  km  over  72  hours. 
Since  the  grid  resolution  is  25  km,  a  3  day  sampling  interval  is  the  minimum 
needed  to  show  expected  movement  of  a  feature  within  the  data  set. 

The  imagery  was  displayed  on  a  Macintosh  II  desktop  computer  using 
Imagic  image  processing  software  (Norris  and  Emery,  1989).  Imagic  was  designed 
to  display  satellite  remote  sensing  imagery  on  the  Macintosh  and  is  distributed  by 
NSIDC  with  the  CD-ROM  data  set.  Images  were  color  enhanced  to  accentuate 
multiyear  ice  features  and  looped  to  provide  a  continuous  animation  of  feature 
movement. 

2.     Feature  Tracking 

The  two  main  steps  for  determining  ice  motion  from  successive  images 
are  tie  point  selection  and  registration.  These  steps  are  achieved  in  this  study 
through  subjective  imagery  analysis  following  rules  that  could  be  incorporated 
into  an  objective  scheme.  Tie  points  are  selected  by  identifying  spatially  coherent 
features  through  a  visual  inspection  of  the  animated  imagery.  These  tie  points  are 
registered  in  succeeding  images  by  matching  the  same  points  in  their  new 
locations.  This  procedure  is  described  in  detail  below: 

1)     Create  a  series  of  multiyear  ice  fraction  data  sets. 
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2)  Load  the  data  sets  into  Imagic  and  generate  ice  fraction  images. 
Decide  on  a  color  enhancement  scheme  that  best  shows  multiyear  ice  features 
and  apply  this  enhancement  to  all  images. 

3)  Using  Imagic,  observe  the  advection  of  the  multiyear  ice  features. 

4)  Select  features  that  appear  to  be  advected  with  the  general  ice  flow 
according  to  the  following  criteria: 

•  Size  -  Features  must  be  small,  no  more  than  100  km  across  (4 
pixels.  Accurate  tiepoint  registration  within  larger  features  is  too 
difficult. 

•  Background  Discrimination  -  Features  must  stand  out  from  the 
background  by  at  least  10  percentage  points  for  most  of  their 
lifespan.  This  increases  one's  confidence  that  what  is  being 
observed  is  actually  a  feature  advected  in  the  flow  rather  than 
random  variability  of  multiyear  ice  fractions. 

•  Temporal  Coherence  -  Features  must  be  present  in  a  least  6 
consecutive  images  (18  days).  This  further  reduces  the  possibility 
of  spurious  tracks  due  to  temporary  brightness  temperature 
fluctuations  in  the  ice  cover. 

5)  Record  the  pixel  position  of  an  appropriate  tie  point  within  the 
feature  during  its  lifespan  using  the  following  guidelines: 

•  If  one  pixel  value  is  3  percentage  points  or  more  higher  than  all 
other  values  within  the  feature,  select  it. 
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•  If  2  or  more  pixel  values  are  within  3  percentage  points,  select 
the  pixel  which  best  corresponds  to  the  previous  relative  position 
of  the  tie  point  within  the  feature. 

•  If  the  relative  maximum  within  the  feature  is  ambiguous, 
choose  the  pixel  which  most  conserves  ice  motion. 

6)  Convert  the  pixel  location  within  the  polar  stereographic  grid  to 
geodetic  latitude  and  longitude  using  software  supplied  by  NSIDC.  A  version 
modified  for  the  Macintosh  is  listed  in  Appendix  A. 

7)  Determine  ice  velocities  by  measuring  the  change  of  position  with 
time.  FORTRAN  code  that  converts  geodetic  positions  to  u  and  v  components  of 
velocity  is  also  contained  in  Appendix  A. 

An  example  of  the  feature  tracking  technique  is  shown  in  Figures  24 
through  26.  This  particular  feature  was  tracked  from  1  February  through  22  April 
1988  as  it  moved  southward  down  the  east  coast  of  Greenland.  Figure  24 
contains  six  consecutive  images  of  multiyear  ice  fraction  contour  plots  from  3-17 
March.  Each  plot  covers  the  same  500  by  250  km  region  centered  at  76.1°N, 
13.5°W  (note  the  Greenland  coast  in  the  lower  left  corner  of  the  plots).  An  area 
of  relatively  high  multiyear  ice  fractions  can  be  seen  migrating  from  the  upper  left 
to  the  lower  right  corner  of  each  plot.  This  is  the  trackable  feature.  Figure  25 
shows  pixel  values  for  this  same  region  expressed  as  multiyear  ice  fractions.  The 
circled  number  is  the  feature  tie  point.  Pixel  locations  were  converted  to  geodetic 
coordinates  and  the  u  and  v  velocity  components  are  shown  in  Figure  26. 
3.     Feature  Selection 

The  number  of  trackable  features  depends  on  the  region  of  interest.  For 
example,  the  generally  homogeneous  conditions  in  the  high  arctic  (above  85°N) 
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create  a  poor  environment  for  feature  tracking.  No  suitable  features  were 
observed  in  the  1988  data  set.  Only  two  features  in  the  high  arctic  were 
observed  during  the  1989  winter  season.  On  the  other  hand,  Fram  Strait  is  a 
potentially  rich  area  due  to  the  merging  of  the  Siberian  and  Polar  branches  of  the 
TDS. 

To  assess  the  validity  of  this  feature  tracking  technique,  the  trajectories 
of  multiyear  ice  features  were  compared  with  true  ice  motion  using  the  arctic  data 
buoy  program  as  ground  truth  for  the  ice  velocities.  The  criterion  to  colocate 
trackable  features  with  a  drifting  buoy  was  challenging.  Suitable  features  had  to 
be  identified  and  their  positions  compared  against  buoy  tracks  as  reported  in 
Colony  and  Rigor  (1990  and  1991).  Three  feature-buoy  matchups  were  selected 
from  the  winter  88  data  set,  each  representing  a  different  region  of  the  arctic. 
Named  after  the  closest  geographic  reference  point,  the  features  are  Severnaya 
Zemyla  (SZ),  Franz  Josef  Land  (FJL),  and  East  Siberian  Sea  (ESS).  A  feature  from 
the  Winter  89  data  set  (CA)  was  chosen  to  represent  motion  in  the  central  arctic 
basin.  Although  no  satisfactory  matchup  between  a  trackable  feature  and  a 
drifting  buoy  was  found  in  the  Fram  Strait  region,  a  buoy  (#7001)  that  passed 
down  the  Greenland  coast  two  weeks  ahead  of  a  feature  (FRAM)  was  used  for  a 
relative  comparison  of  the  ice  trajectory. 

The  tracked  features  and  ground  truth  buoys  were  not  coincident  in 
space;  the  average  distance  between  the  selected  features  and  buoys  varied  from 
76  to  557  km.  The  closest  matchup  (76  km)  occurred  between  FJL  and  buoy 
#1790.  The  ESS  feature  was  420  km  from  buoy  #1896.  SZ  was  320  km  from  the 
buoy  #7048.  CA  averaged  497  km  from  the  closest  buoy  (#3288).  Table  4 
summarizes  the  characteristics  of  the  feature-buoy  tracks. 
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TABLE  4:  MYI  FEATURE  -  BUOY  DATA 


Feature 

Buoy 

Separation 
(km) 

Start  Track 
Date 

Position 

Stop  Track 
Date 

(Total  Days) 

Position 

ESS 

7048 

420 

1  APR  88 

76.2N, 
160.1E 

25  APR  88 
(24) 

77.3N, 
159.0E 

FJL 

1790 

76 

1  FEB  88 

80.2N, 
41.3E 

14  MAR  88 
(42) 

79.5N, 
41.9E 

SZ 

1897 

320 

1  FEB  88 

81.7N, 
105.9E 

28  APR  88 
(85) 

82.0N, 
100.6E 

CA 

3288 

491 

1  JAN  89 

84.5N, 
143.5E 

31  MAR  89 
(90) 

86.6N, 
125.2E 

FRAM 

N/A 

N/A 

1  FEB  88 

79.9N, 
4.6W 

22  APR  88 
(80) 

68.1N, 
23 .9E 

Geolocation  Key:   SZ  -  Sevemaya  Zemlya 

FJL  -  Franz  Josef  Land 

ESS  -  East  Siberia  Sea 

CA  -  Central  Arctic 

FRAM  -Fram  Strait 


Note:  Features  are  named 
according  to  nearby 
geographic  reference 
points. 


4.     Errors 

Before  ice  motion  vectors  between  the  buoy  -  feature  pairs  are 
compared,  the  maximum  errors  inherent  to  this  technique  are  described.  Three 
types  of  positional  errors  are  possible:  imperfect  tiepoint  registration,  SSM/I  data 
geo-location,  and  ARGOS  buoy  positioning.  In  addition,  ice  velocity  gradients 
between  the  feature  and  buoy  locations  need  to  be  considered.  These  errors, 
expressed  as  an  ice  speed  error  over  a  3-day  period,  are  examined  in  detail  below. 

Imperfect  tiepoint  registration  within  the  data  grid  is  the  largest  source 
of  error  in  this  feature  tracking  technique.  For  small,  well-defined  features  the 
maximum  registration  error  is  estimated  to  be  one  pixel  for  each  observation.  The 
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actual  error  should  be  uniformly  distributed  within  the  pixel.  The  maximum  error 
is  equivalent  to  a  50  km  or  19.3  cm/s  velocity  error  over  a  3  day  period.  The  geo- 
location  error,  as  reported  by  NSIDC  (1990),  is  8  km  for  each  scan.  Since  each 
pixel  represents  a  daily  averaged  mean  of  multiple  scans,  the  geo-location  error 
should  be  somewhat  smaller.  For  an  error  of  5  km,  the  magnitude  of  ice  motion 
could  vary  by  3.8  cm/s  over  3  days.  Thus  the  maximum  estimated  position  error  is 
30  km.  If  this  error  is  assumed  to  be  uniformly  distributed,  then  the  standard 
deviation  of  the  error  is  19.0  km. 

The  estimated  velocity  error  due  to  the  maximum  position  error  is  23.1 
cm/s  over  3  days.  This  is  over  twice  the  typical  ice  velocity  of  10  cm/s.  If  instead 
the  standard  deviation  of  the  position  error  is  considered,  the  velocity  error 
reduces  to  14.7  cm/s.  These  velocity  errors  decrease  if  the  motion  is  averaged 
over  longer  time  periods.  For  example,  a  12  day  average  would  result  in  a 
maximum  error  of  5.8  cm/s. 

The  ARGOS  system  currently  used  to  track  buoy  locations  can 
determine  position  accuracy  within  a  few  hundred  meters  (Thorndike,  1986). 
Assuming  a  maximum  error  of  500  m  for  each  ARGOS  reported  position,  the 
resulting  velocity  error  over  3  days  would  be  0.4  cm/s. 

When  comparing  SSM/I-derived  ice  trajectories  with  buoy  locations,  the 
velocity  gradient  between  the  two  locations  must  be  considered.  These  velocity 
gradients  are  the  result  of  the  ice  pack  response  to  external  forces  such  as  wind 
and  ocean  currents.  Thorndike  (1986)  compared  the  velocity  components  of  ice 
camps  and  buoys  separated  by  distances  ranging  from  1  to  1000  km.  He 
concluded  that  on  scales  of  this  magnitude  ice  behaves  as  a  homogeneous, 
isotropic  velocity  field.    Thorndike  examined  the  correlations  of  buoy  and  ice 
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camp  velocities  as  a  function  of  distance.  For  a  100  km  camp-buoy  separation,  he 
found  the  correlation  to  be  0.9.  At  a  400  km  separation,  the  correlation  drops  to 
approximately  0.6. 

Using  these  correlations,  the  variance  of  the  differences  in  ice  velocity 
between  two  separate  points  within  the  ice  pack  can  be  estimated.  The  variance 
of  the  velocity  difference  between  the  points  (a2err)  is  given  by 

a2err  =  E[V(x+s)-V(x)]2  (8) 

where  E  is  the  expected  value  operator  and  V(x)  and  V(x-i-s)  are  the  ice  velocities 
separated  by  some  distance  s.  Assuming  the  variance  of  ice  velocity  is  the  same 
at  both  points,  (  8)  reduces  to 

a2err  =  2o2i(i-p)  (9) 

where  G2i  is  the  variance  of  the  ice  velocity  and  p  is  the  correlation  coefficient. 
Given  a  standard  deviation  of  ice  velocity  of  10  cm/s  (Thorndike,  1986)  and 
assuming  the  u  and  v  components  of  motion  are  independent  of  each  other  (the 
total  variance  then  equals  the  sum  of  the  variance  of  each  component),  <72err  is  40 

cm2/s2  and  10  cm2/s2  for  correlations  of  0.6  and  0.9,  respectively.  The  resulting 
velocity  errors  are  6.3  cm/s  for  a  400  km  separation  and  3.2  cm/s  for  a  100  km 
separation.  These  velocity  errors  equate  to  16  km  and  8  km  positional  errors, 
respectively. 

Table  5  summarizes  the  maximum  errors  expected  for  each  error  source. 
The  first  three  position  errors  are  per  observation.    The  corresponding  velocity 
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error  results  from  combining  two  position  errors  for  each  velocity  calculation. 
Velocity  errors  are  calculated  over  3  and  12  day  day  intervals  to  show  how  the 
maximum  estimated  error  decreases  for  longer  mean  flow  calculations. 


TABLE  5:  SUMMARY  OF  MAXIMUM  ERRORS 


Error  Source 

Positional  Error 
(km) 

3 -Day 

Velocity  Error 

(cm/s) 

12  -  Day 

Velocity  Error 

(cm/s) 

Tie  Point 
Registration 

25 

19.3 

4.8 

SSM/I 
Geolocation 

5 

3.8 

0.9 

Buoy 
Position 

0.5 

0.4 

0.1 

Velocity 
Gradient 
(100  Km) 

8  (3  Day) 

3.2 

3.2 

Velocity 
Gradient 
(400  Km) 

16  (3  Day) 

6.3 

6.3 

This  error  analysis  shows  that  3  day  ice  motion  estimates  are  unreliable 
indicators  of  true  ice  motion  because  the  velocity  error  can  exceed  the  typical 
true  ice  velocities  of  10  cm/s.  Longer  term  averages,  however,  reduce  the  impact 
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of  the  SSM/I  position  errors  and  have  the  potential  to  represent  accurately 
general  ice  flow. 

An  expected  correlation  coefficient  can  be  calculated  from  the  maximum 
expected  errors.  This  correlation  can  be  used  to  validate  the  hypothesis  that  the 
feature  tracks  are  representative  of  true  ice  motion.  If  the  observed  correlations 
are  less  than  the  expected  value,  then  either  the  feature  tracks  do  not  represent 
true  ice  motion  or  the  errors  were  underestimated.  If  the  observed  correlations  are 
larger  than  the  expected  value,  then  the  actual  errors  are  smaller  than  expected. 

The  correlation  coefficient  between  two  quantities,  x  and  y,  is  given  by 

Exy  /im 

P  =  7^ — F  \\n  (10) 

[Exx  EyyJ1/z- 

where  both  x  and  y  are  assumed  to  have  zero  means.  Let  x  and  y  be  defined  as 
the  ice  motion  at  the  feature  and  buoy,  respectively,  as  observed  from  the  data: 

x  =  v  +  ef+ed  (ll) 

y  =  v  +  £b  (12) 

where  v  is  the  true  ice  velocity,  ef  is  the  error  in  measuring  the  velocity  of  the 
feature,  ed  is  the  error  due  to  the  velocity  difference  between  the  feature  and  the 
buoy  and  eb  is  the  error  in  measuring  the  velocity  of  the  buoy.  The  variances  due 
to  these  errors  are  defined  as 

o2f  =  Eef2  +  E£d2  (13) 

o2b  =  Eeb2  (14) 
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a2i  =  Evz  (15) 

where  a-f  is  the  variance  caused  by  the  feature  position  and  velocity  difference 
errors,  ozb  is  the  variance  due  to  the  buoy  position  errors  and  Gzi  is  the  variance 
of  true  ice  velocity.  The  expected  values  in  equation  (10)  become 


Exy  =  ozi  (16) 

EXx  =  azi  +  azf  (17) 

Eyy  =  azi  +  ozb  (18) 


Thus  the  expected  correlation  coefficient  can  be  defined  as  the  variance 
of  true  ice  motion  divided  by  the  square  root  of  the  variances  of  buoy  motion  and 
feature  motion  multiplied  together: 

°2i 

(19) 


[(oVo2b)(a2i  +  a2f)]l/2 


As  before,  azj  for  each  component  of  motion  is  50  cmz/sz.  Using  the  3  day 
sampling  interval  from  Table  5,  c\  is  negligible  at  0.16  cmz/sz  and  ozf  is  the  sum 
of  the  variances  due  to  positional  error  (266.8  cmz/s2)  and  the  velocity  gradient 
between  the  buoy  and  the  feature  (39.7  cmz/sz  for  400  km  and  10.2  cm2/sz  for 
100  km  separations).  Using  these  values  in  equation  (19),  the  expected 
correlation  for  a  feature  and  buoy  100  km  apart  is  0.391.  For  a  separation  of  400 
km,  the  correlation  decreases  to  0.374. 
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These  numbers  indicate  that  if  the  correlation  between  3  day  mean 
velocities  determined  from  buoys  and  multiyear  ice  features  exceeds 
approximately  0.4,  the  feature  is  moving  with  the  ice  and  the  errors  associated 
with  the  tracking  technique  were  less  than  expected.  If  the  same  calculations  are 
performed  for  12  day  mean  ice  velocities,  correlations  rise  to  0.866  and  0.783  for 
100  and  400  km  separations,  respectively. 

C.   RESULTS 

1.      Ice  Trajectories 

Trajectories  of  the  four  feature  -  buoy  pairs  from  the  different  arctic 
regions,  each  plotted  at  3  day  intervals,  are  shown  in  Figures  27  through  31. 
Because  the  feature  tracks  are  limited  to  discrete  25  km  jumps  in  position,  periods 
of  slow,  looping  or  retrograde  motion  confuse  the  trajectories.  In  these  cases,  the 
positions  have  been  slightly  offset  to  show  movement  more  clearly. 

The  correspondence  between  3  day  positions  of  the  buoys  and  features 
appears  somewhat  erratic.  This  is  probably  due  to  the  resolution  limitations 
inherent  in  the  SSM/I  data  grid.  Taken  over  longer  periods,  however,  the  general 
flow  corresponds  quite  well  to  the  buoy  positions.  This  is  best  shown  in  the  Fram 
Strait  trajectories  (Figure  27).  Both  the  buoy  and  FRAM  tracks  are  roughly 
parallel  as  they  traveled  southward  down  the  east  Greenland  coast,  the  buoy 
preceding  FRAM  by  about  two  weeks.  Of  special  note  is  the  sudden  surge 
eastward  at  78°N  for  the  feature  and  at  69°N  for  the  buoy.  Although  separated 
by  9°  of  latitude,  this  surge  occurred  at  the  same  time  (20-23  Feb).  This 
corresponds  to  a  period  of  strong  geostrophic  wind  flow  to  the  northeast  due  to  a 
984  mb  low  centered  over  the  northern  coast  of  Greenland  (Colony  and  Rigor, 
1990). 
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2.      Velocity  Comparisons 

For  each  feature  -  buoy  pair,  the  u  and  v  components  of  motion  over 
each  3  day  period  were  calculated.  Buoy  velocities  were  computed  directly  from 
their  reported  positions.  With  the  exception  of  ESS,  all  multiyear  feature 
velocities  were  smoothed  over  two  positions  to  reduce  velocity  spikes  resulting 
from  the  gridded  data  set.  The  ESS  feature,  tracked  for  24  days,  has  only  one 
local  velocity  maximum  in  each  direction  so  the  smoothing  was  not  applied. 

Correlation  coefficients  between  the  buoys  and  both  smoothed  and 
unsmoothed  feature  velocities  are  presented  in  Table  6.  Applying  a  t  distribution 
test  to  the  smoothed  velocities,  all  but  one  correlation  (Franz  Josef  Land  u 
component)  exceed  the  95%  confidence  level  that  the  feature  and  buoy 
velocities  were  correlated.  The  generally  high  correlations  are  particularly 
surprising  since  the  maximum  error  analysis  indicates  that  correlations  only  as 
high  as  0.4  should  be  expected  for  3  day  mean  ice  velocities.  These  higher  than 
expected  correlations  imply  that  the  estimated  errors  used  in  the  error  analysis 
were  too  large. 

Specific  velocity  comparisons  are  plotted  in  Figures  32  through  35  and 
are  discussed  below. 

a.    Franz  Josef  Land  -  Buoy  1790  (Figure  32) 

This  is  the  closest  spatial  match-up  of  any  feature  -  buoy  pair 
investigated.  The  separation  distance  averaged  76  km  over  a  42  day  period.  The 
buoy  trajectory  shows  complex  motion  with  two  separate  periods  of  retrograde 
east  -  west  movements.  The  SSM/I  grid  had  trouble  resolving  these  complicated 
movements,  providing  the  lowest  correlation  of  any  velocity  component:  -0.069 
for  the  unsmoothed  and  0.322  for  the  smoothed  u  component.    When  velocities 
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were  averaged  over  12  days,  the  unsmoothed  u  component  correlation  improved 
to  0.666.  Three  day  v  correlations  were  0.458  and  0.678  for  unsmoothed  and 
smoothed  velocities,  respectively.  The  higher  correlations  in  the  v  direction  are 
probably  due  to  the  slow,  southward  trend  of  the  drift. 


TABLE  6:  MYI  FEATURE  -  BUOY  CORRELATIONS 

Feature 

Buoy 

Correlations  (3  day  velocities) 

u 

V 

smooth  u 

smooth  v 

ESS 

1896 

.791 

.724 

.845 

.490 

FJL 

1790 

-.069 

.458 

.322 

.674 

SZ 

7048 

.652 

.537 

.685 

.455 

CA 

3288 

.051 

.501 

.311 

.582 

b.    East  Siberian  Sea  -  Buoy  1896  (Figure  33) 

This  pair  was  tracked  over  the  shortest  interval  (24  days).  Although 
separated  by  320  km,  it  also  had  the  highest  correlations.  The  correlations  for  the 
unsmoothed  u  and  v  components  were  0.791  and  0.724,  respectively.  The 
trajectory  characteristics  responsible  for  these  high  correlations  were  a  steady 
northwestern  drift  interrupted  by  a  strong  southwesterly  surge  (16-19  April)  large 
enough  to  be  resolved  within  the  SSM/I  grid.  This  surge  corresponded  to  a 
period  of  strong  geostrophic  wind  flow  to  the  southwest  from  a  low  in  the 
Chukchi  Sea  region  (Colony  and  Rigor,  1990). 
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c.  Severnaya  Zemyla  -  Buoy  7048  (Figure  34) 

The  feature  and  buoy  were  tracked  for  75  days.  The  average 
separation  distance  was  420  km.  Both  trajectories  showed  a  slow  easterly  drift 
culminating  in  a  clockwise  loop.  Correlations  for  the  smoothed  u  and  v 
components  were  0.685  and  0.455,  respectively.  The  higher  u  correlation  is  due 
to  the  predominantly  eastward  movement  of  the  feature  and  buoy. 

d.  Central  Arctic  -  Buoy  3288  (Figure  35) 

This  pair  was  taken  from  the  1989  data  set  to  represent  ice  motion  in 
the  central  arctic  region.  They  were  tracked  for  90  days  at  an  average  separation 
of  491  km.  The  smoothed  u  correlation  was  low  at  0.31 1.  The  v  component  was 
much  higher  at  0.582.  The  low  u  correlation  can  be  explained  by  the  steady 
northerly  movement  of  the  feature.  East-west  velocities  were  generally  small  and 
not  well  resolved  in  the  SSM/I  data  grid.  Correlations  for  the  12  day  mean 
velocities  jumped  to  0.783  and  0.776  for  u  and  v  components,  respectively.  The 
longer  sampling  period  allowed  time  for  the  slow  motion  of  the  ice  field  in  this 
region  to  be  resolved  within  the  SSM/I  grid. 
3.      Position  Error 

The  maximum  estimated  feature  position  error  can  be  compared  against 
the  position  errors  calculated  from  the  data.  Solving  equation  (19)  for  a-f: 


a2f=a2j77-ll  (20) 


.9        _?  (  * 


a^i  is  assumed  to  be  the  same  as  the  variance  of  buoy  velocities  and  a-b  is 
negligible  compared  with  o^ice.    The  average  variance  of  each  buoy  velocity 
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component  is  30.08  cm^/s2  and  the  average  correlation  (p)  between  the 
smoothed  feature  and  buoy  velocities  is  0.545. 

Using  equation  (20)  and  the  average  variances  and  correlations 
observed  from  the  feature  track  data,  o2f  is  71.19  cm2/s2.  The  ice  velocity 
difference  between  the  buoy  and  feature  must  be  taken  into  account  because  c2f 
is  the  sum  of  the  variances  due  to  position  error  and  the  velocity  difference. 
Substituting  the  average  standard  deviation  of  the  buoy  velocities  into  equation 
(9)  and  again  assuming  that  u  and  v  components  of  motion  are  independent  of 
each  other,  the  variance  due  to  the  ice  velocity  gradient  is  12.32  cm2/s2  and  3.08 
cm2/s-  for  400  km  and  100  km  feature-buoy  separations,  respectively.  Taking 
into  account  the  average  feature-buoy  separation  was  326  km,  nearly  17%  of  the 
observed  variance  can  be  explained  by  the  ice  velocity  gradient  between  the 
feature  and  the  buoy.  The  remaining  variance  of  approximately  58  cm2/s2 
equates  to  a  10.7  cm/s  total  velocity  error  (u  and  v  velocity  variances  added 
together). 

This  velocity  error  is  less  than  half  of  the  maximum  expected  value  of 
23.1  cm/s  and  72%  of  the  standard  deviation  of  the  expected  error.  A  10.7  cm/s 
velocity  error  over  a  three  day  period  implies  close  to  a  15  km  position  error  in 
each  of  the  two  observations  used  to  calculate  the  velocity.  This  is  less  than  both 
the  maximum  expected  postion  error  (30  km)  and  the  standard  deviation  of  this 
error  (19  km).  Possible  explanations  for  the  difference  between  expected  and 
actual  values  include:  1)  the  tiepoint  registration  and  geo-location  errors  are 
smaller  than  expected  because  the  gridded  brightness  temperatures  in  the  dataset 
are  obtained  by  averaging  information  gathered  from  several  satellite  passes  each 
day;  2)  the  actual  ice  velocity  gradient  between  the  features  and  buoys  is  smaller 
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than  Thorndike's  (1986)  estimates;  and  3)  the  standard  deviation  of  the  ice 
velocities  is  smaller  than  the  assumed  value  of  10  cm/s.  The  latter  explanation 
probably  had  the  greatest  influence  because  standard  deviations  of  the  buoy 
velocities  only  ranged  from  3  to  9  cm/s.  These  smaller  than  expected  standard 
deviations  could  result  from  the  the  use  of  three  day  ice  velocity  averages  which 
reduce  the  effect  of  short  term  velocity  surges. 

D.   DISCUSSION 

There  is  good  reason  to  doubt  the  practicality  of  inferring  sea  ice  motion  by 
tracking  multiyear  ice  concentrations  from  space-based  passive  microwave 
radiometers.  Data  resolution  is  poor.  At  37  and  19.35  GHz  (the  frequencies 
currently  used  to  calculate  ice  fractions),  resolutions  are  33  and  56  km, 
respectively.  The  25  km  grid  resolution  of  the  SSM/I  dataset  is  only  achieved  by 
averaging  several  observations  over  the  same  spot  each  day.  Not  only  does  the 
low  resolution  increase  the  position  error,  it  also  greatly  reduces  the  possibility  of 
tracking  individual  floes.  The  radiometer  must  instead  observe  broad  regions  of 
anomalous  radiances  to  create  trackable  features.  The  precision  of  the  multiyear 
ice  fraction  algorithm  is  only  estimated  to  be  in  the  range  of  13  -  25  %.  Moreover, 
a  variation  of  0.04  cm  in  cloud  liquid  water  content  will  contribute  as  much  as 
30%  to  the  uncertainty  of  computed  multiyear  ice  fractions  (Cavalieri  et  al., 
1984). 

In  spite  of  these  handicaps,  ice  drift  inferred  through  tracking  multiyear  ice 
features  from  SSM/I  radiometer  data  corresponded  quite  well  with  ice  motion  as 
measured  by  in-situ  buoys.  Ice  feature  trajectories  were  similar  to  nearby  buoy 
drift  patterns.  In  two  cases  (ESS  and  FRAM),  when  high  atmospheric  pressure 
gradients  caused  anomalously  strong  wind  stress  on  the  ice  pack,  both  the  feature 
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and  buoy  exhibited  similar  trajectories  in  accordance  with  the  forcing. 
Comparisons  between  the  orthogonal  velocity  components  of  the  feature  and 
buoy  motion  averaged  over  3  days  also  showed  surprisingly  high  correlations. 
Instances  of  high  correlations  in  one  component  direction  but  low  correlations  in 
the  other  can  be  explained  by  slow  or  oscillatory  movement  in  the  direction  of 
the  velocity  component  exhibiting  the  low  correlation  and  are  a  function  of  the 
low  spatial  resolution  of  the  grid. 

The  good  results  achieved  above  were  the  result  of  being  extremely  selective 
during  the  tiepoint  selection  process.  First,  all  tracked  features  were  of  high 
multiyear  ice  fractions  surrounded  by  a  background  of  relatively  low  fractions. 
Multiyear  ice  is  thicker  than  first  year  ice  and  therefore  less  subject  to 
deformational  and  destructive  processes  that  might  obscure  the  characteristic 
radiometric  signal.  Second,  the  features  were  small — no  more  than  100  km  across. 
This  simplified  the  task  of  registering  the  same  tiepoint  in  successive  images.  The 
variation  of  multiyear  ice  fraction  percentages  within  a  larger  feature  could 
preclude  accurate  tiepoint  registration.  Limiting  the  feature  spatial  extent  also 
reduced  the  possibility  of  confusing  an  atmospheric  disturbance  with  an  ice 
feature.  Third,  imposing  a  temporal  coherence  requirement  of  18  days  ensured 
that  the  tracked  object  was  not  a  result  of  random  brightness  temperature 
fluctuations  in  the  ice  cover.  Finally,  limiting  the  record  to  the  winter  season 
reduced  the  possibility  of  surface  melt  conditions  contaminating  radiances  from 
the  sea  ice.  Another  benefit  of  the  winter  season  is  that  moisture-laden  clouds 
with  their  contaminating  effects  are  generally  restricted  to  the  open  ocean  and 
the  marginal  ice  zone  areas. 


66 


The  number  of  trackable  features  was  limited.  The  abovementioned 
limitations  restricted  the  number  of  suitable  features  to  only  a  dozen  or  so  during 
the  winter  season.  Relaxing  these  restrictions  should  produce  more  trackable 
features  (and  possibly  introduce  spurious  features)  and  is  an  area  suitable  for 
further  research.  Development  of  an  objective  scheme  to  identify  and  track' 
features  should  also  increase  the  number  of  available  datapoints. 

Despite  the  above  restrictions,  this  limited  ability  to  observe  and  track  areas 
of  anomalous  multiyear  ice  fractions  can  be  useful  for  arctic  studies.  Ice  fluxes 
from  the  central  arctic  have  been  shown  to  be  the  largest  exporter  of  fresh  water 
from  the  Arctic  Ocean  and  the  largest  freshwater  input  to  the  Greenland  Sea 
(Aagaard  and  Carmack,  1989).  Available  databases,  extending  back  to  1978,  can 
be  used  to  assess  the  intra  and  interannual  variability  of  ice  transport. 
Researchers  can  identify  previous,  and  possibly  forecast  future,  pulses  of  ice 
types  to  examine  their  influence  on  regional  heat  and  salinity  budgets.  Ice  camps 
are  generally  located  on  relatively  smooth,  thick  floes  of  multiyear  ice  for 
logistical  and  safety  reasons.  Potential  sites  are  currently  scouted  by  aircraft. 
This  space-based  technique  offers  the  possibility  of  remotely  identifying  sites  well 
in  advance  of  deployment.  As  satellite  instrumentation  improves  and  emanations 
from  higher  microwave  frequencies  are  utilized,  surface  resolution  will  also 
improve  resulting  in  increased  reliability  of  short  term  ice  velocity  vectors. 

E.    CONCLUSIONS 

When  images  of  arctic  multiyear  ice  fractions  calculated  from  winter  season 
SSM/I  data  are  examined  in  sequence,  distinct  features  are  observed  to  move  in 
the  general  direction  of  climatological  ice  flow.  These  features  exist  throughout 
the  arctic  basin  but  are  mainly  located  in  regions  where  first  year  and  multiyear 
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ice  merge  together.  Some  features  are  short-lived,  others  can  be  tracked 
throughout  the  data  record. 

Four  of  these  multiyear  ice  features  from  the  winters  of  1988  and  1989  were 
compared  against  nearby  arctic  buoy  trajectories.  Three  day  average  ice  velocity 
vectors  were  found  to  correlate  well  with  the  buoy  velocities;  the  average 
smoothed  velocity  correlation  was  0.545.  The  maximum  estimated  velocity  errors, 
mainly  due  to  the  low  data  resolution,  are  double  a  typical  ice  velocity  of  10  cm/s 
making  short-term  ice  motion  calculations  unreliable.  Actual  velocity  errors 
derived  from  the  feature-buoy  correlations  are  half  the  maximum  estimated  error 
and  less  than  three  quarters  of  the  standard  deviation  of  this  error.  This  difference 
is  probably  due  to  the  smaller  than  estimated  variances  of  ice  motion  because  the 
ice  motion  is  averaged  over  a  three  day  period,  smoothing  out  short-term  ice 
velocity  surges. 

The  higher  than  expected  correlations  between  multiyear  ice  feature 
velocities  and  arctic  buoys  are  due  to  the  stringent  requirements  established  for 
feature  selection.  Relaxing  these  requirements  could  substantially  increase  the 
number  of  trackable  features  from  the  dozen  or  so  noted  in  the  winter  of  1988  but 
also  decrease  the  confidence  that  tracked  features  accurately  represent  ice 
motion.  A  similar  increase  in  trackable  features  could  probably  be  attained  by 
objectively  scanning  the  dataset  for  subtle  features  that  were  not  recognized 
through  a  subjective  analysis  of  the  data. 
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Figure  22       Field  of  mean  ice  motion  (from  Colony  and  Thorndike,  1984). 
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Figure  23       Average   sea  ice   drift  showing   the   two  components   of  the 
Transpolar  Drift  Stream  (  from  Polarstern  Scientific  Party,  1988). 
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Figure  24a  Contour  plots  showing  advection  of  multiyear  ice  feature  along  the 
east  Greenland  coast.  Images  from  2  and  5  March,  1988.  Plots  are 
centered  at  76.1°N,  13.5°W.  The  feature  moves  from  the  upper  left 
to  the  lower  right  corner. 
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Figure  24b     Same  as  Figure  24a  except  images  are  from  8  and  1 1  March,  1988. 
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Figure  24c      Same  as  Figure  24a  except  images  are  from  14  and  17  March,  1988. 
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Figure  25a     Pixel  values  for  multiyear  ice  contour  plot  shown  in  Figure  24a.    Circli 
value  is  tiepoint  registration  location. 
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Figure  25b 


Pixel  values  for  multiyear  ice  contour  plot  shown  in  Figure  24b.    Circled 
value  is  tiepoint  registration  location. 
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Figure  25c      Pixel  values  for  multiyear  ice  contour  plot  shown  in  Figure  24c. 
Circled  value  is  tiepoint  registration  location. 
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Figure  26       Ice  velocities  for  Fram  Strait  feature  (1  February  -  22  April,  1988). 
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Figure  27       Ice  trajectories:  Fram  Strait  feature  and  Buoy  7001  (1  February  -  22 
April,  1988).  Circled  points  are  initial  positions. 
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Figure  28       Ice  trajectories:     Franz  Josef  Land  feature  and  Buoy   1790  ( 
February.  -  14  March,  1988).  Circled  points  are  initial  positions. 
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Figure  29       Ice  trajectories:    East  Siberia  Sea  feature  and  Buoy  1896  (1 
April,  1988).  Circled  points  are  initial  positions. 
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Figure  30       Ice  trajectories:     Severnaya  Zemyla  feature  and  Buoy  7048  ( 
February  -  28  April,  1988).  Circled  points  are  initial  positions. 
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Figure  31        Ice  trajectories:   Central  Arctic  feature  and  Buoy  3288  (1  January 
31  March,  1988).  Circled  points  are  initial  positions. 
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Figure  32       Franz  Josef  Land  and  Buoy  1790  velocity  comparison.    Buoy  is 
solid  line,  feature  is  dashed  line. 
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Figure  33       East  Siberian  Sea  and  Buoy  1896  velocity  comparison.    Buoy  is 
solid  line,  feature  is  dashed  line. 
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Figure  34       Severnaya  Zemyla  and  Buoy  7048  velocity  comparison.    Buoy  is 
solid  line,  feature  is  dashed  line. 
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Figure  35       Central  Arctic  and  Buoy  3288  velocity  comparison.    Buoy  is  solid 
line,  feature  is  dashed  line. 
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V.   CONCLUSIONS  AND  RECOMMENDATIONS 

Arctic  sea  ice  conditions  are  an  important  component  of  the  global  climate 
system.  The  sea  ice  cover  affects  the  Earth's  radiative  balance  by  influencing  the 
heat  flux  between  the  ocean  and  atmosphere  and  heat  transport  between  the 
lower  and  upper  latitudes.  Passive  microwave  remote  sensing  of  sea  ice  offers 
continuous,  all-weather  observations  of  these  remote  areas  independent  of  solar 
illumination.  This  is  particularly  important  in  the  arctic  where  much  of  the  region 
is  shrouded  by  clouds  or  covered  in  darkness. 

Accurate  ice  observations  using  passive  microwave  remote  sensing 
techniques  require  knowledge  of  the  radiative  characteristics  of  sea  ice.  The 
radiative  properties  of  sea  ice  are  complicated  due  to  its  complex  nature  which  is 
determined  by  formation  environment,  age  and  season.  Other  complicating 
factors  include  contamination  from  moisture-laden  clouds  and  the  presence  of 
snow  cover,  melt  ponds  and  ridging.  In  spite  of  these  difficulties,  emissivity 
differences  between  sea  ice  and  water  at  various  microwave  frequencies  can  be 
exploited  to  estimate  ice  characteristics.  Current  algorithms  use  polarization  and 
gradient  ratio  differences  to  estimate  ice  concentrations  and  type. 

Historical  SSM/I  data  are  readily  available  on  CD-ROM  media  distributed  by 
NSIDC.  Two  applications  of  this  dataset  to  arctic  sea  ice  studies  were  examined. 
Data  from  1  February  to  15  October  1988  were  used  to  associate  anomalous 
diurnal  tidal  currents  over  the  Yermak  Plateau  with  ice  morphology.  Data  from  1 
February  to  28  April  1988  and  1  January  to  31  March  1989  were  used  to  explore 
the  possibility  of  determining  sea  ice  motion  by  tracking  advection  of  multiyear 
ice  features.   Software  supplied  by  NSIDC  were  modified  to  compute  and  display 
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the  requisite  ice  concentration,  gradient  ratio  and  multiyear  ice  fraction 
calculations  and  imagery  on  a  Macintosh  II  desktop  computer. 

The  hypothesized  ice  concentration/morphology  changes  induced  by  the 
anomalously  strong  diurnal  tidal  currents  observed  over  the  Yermak  Plateau 
could  not  be  verified  by  the  SSM/I  data.  Either  the  currents  do  not  affect  the  ice 
cover  or  SSM/I  limitations  obscure  the  ice  signatures  of  the  currents.  The  low 
spatial  resolution  and  susceptibility  to  melt  and  snow  contamination  of  passive 
microwave  emissions  prevent  a  definitive  assessment  of  tidal  current  effects  on  ice 
cover  in  this  region. 

Distinct  multiyear  ice  features  are  observed  to  move  in  the  general  direction 
of  the  mean  ice  flow.  Feature  movement  correlates  well  with  ice  motion  as 
measured  by  arctic  data  buoys.  The  subjective  nature  of  the  feature  tracking 
technique  and  its  reliability  on  'features  of  opportunity'  prevent  wide  scale 
application  throughout  the  arctic  region.  Positional  errors  resulting  from  low  data 
resolution  could  make  short  term  ice  motion  calculations  unreliable.  Long  term 
averages,  however,  are  reliable  and  can  be  used  to  augment  our  knowledge  of 
general  ice  flow.  For  example,  motion  vectors  can  be  used  to  estimate  the 
interannual  variability  of  ice  transport  through  the  Fram  Strait  and  aid  our 
understanding  of  heat  and  salinity  budgets  in  the  Greenland  Sea.  In  spite  of  the 
limitations,  tracking  multiyear  ice  features  can  be  a  useful  method  of  discerning 
continuous,  all-weather  ice  motion  from  satellite  data. 

Further  attempts  to  verify  anomalous  diurnal  currents  over  the  Yermak 
Plateau  using  SSM/I  data  are  not  recommended.  The  data  resolution  is  too  low  to 
observe  small  scale  changes  of  sea  ice  concentration  or  morphology.  Feature 
tracking,  however,  does  show  some  promise.     Conversion  of  the  subjective 
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technique  described  here  to  an  objective  scheme  for  computer  analysis  would 
greatly  reduce  the  time  and  effort  required  for  tie  point  selection  and  registration, 
and  probably  increase  the  number  of  tie  points  suitable  for  tracking.  Exploration 
of  the  higher  spatial  resolutions  afforded  by  85.5  GHz  data  would  greatly 
increase  the  reliability  of  ice  motion  vectors. 
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APPENDIX  A 

This  appendix  contains  FORTRAN  source  code  listings  of  software  used  to 
process  the  dataset  contained  on  the  NSIDC  CD-ROM  disks. 


Program 


ssmnce 


locate 


u,v  calc 


Subroutine 


usetie 


ice25 


Page 
89 

95 

102 


Page 
105 

108 


Purpose 


Computes  ice  concentrations  from 
SSM/I  brightness  temperature  values  on 
the  NSIDC  CD-ROM  dataset. 

Transforms  I,  J  coordinates  of  the  NSIDC 
SSM/I  polar  grids  to  latitude  and 
longitudes  coordinates  and  vice  versa. 

Calculates  u  and  v  components  of 
motion  from  a  list  of  user  specified 
geodetic  latitude  and  longitude 
positions. 

Purpose 

Computes  coefficients  for  the  ice 
algorithm  in  ssmiice  using  NASA- 
approved  tiepoints. 

Computes  ice  concentrations  from 
NSIDC  CD-ROM  data  at  25  km  grid 
resolution. 
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program  ssmnce 
cccccccccccccccccccccccccccccccccccccx:ccccccccccx:cccccccccccccccccccc 
c 

c  SSMIICE  -  This  program  computes  ice  concentrations  from  SSMI  TB 
c  grids  residing  on  CD-ROM. 

c 

c     Inputs 
c 
c       USER  DEFINED 

c 

c 

c         drive    :  Name  of  CD-ROM  disk 

c 

c        indate  :  Date  in  the  format  YYMMDD  (eg.,  870710) 

c  From  this  variable  the  CD-ROM  is  searched  for  the 

c  appropriate  filename.  The  filename  for  the  ice 

c  concentration  algorithm  is  derived  from  this 

c  variable. 

c 

c        ihem     :  Defines  the  appropriate  hemisphere  to  apply  this 

c  algorithm  (l=North,  2=South) 

c 

c        icepar  :  Defines  the  type  of  ice  concentrations  to  derive: 

c  l=total,  2=multi-year,  3=first-year. 

c 

c       INTERNAL 

c 

c  infile  :  The  path,  name,  and  descriptor  string  describing 

c  the  data  file  from  which  the  SSMI  Tb  data  will  be 

c  extracted  from  the  CD-ROM.  User  must  specify  appropriate 

c  filepath. 

c 

c  buff    :  Input  buffer  containing  the  SSMI  Tb  data  read  from 

c  the  CD-ROM.  Two  lines  of  SSMI  TB  grid  data  will 

c  read  with  a  single  10  request. 

c 

c  mask     :  Landmask  used  in  computing  the  25  km 

c  ice  concentration  grid. 

c 

c  idiml    :  Number  of  parameters  in  an  SSM/I  Tb  grid  cell 

c 

c  idim2    :  Number  of  cells  in  the  SSM/I  Tb  grid  line 

c 

c  idim3    :  Number  of  grid  lines  in  an  SSM/I  Tb  grid 

c 

c  grdtyp  :  Three  character  descriptor  designating  the  grid  type 

c  to  access  from  the  CD-ROM  (N3B  or  S3B) 
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c        icetyp  :  Three  character  descriptor  designating  the  extension 

c  for  the  resultant  ice  parameter  image  (N3C  or  S3C). 

c 

c        c        :  Array  of  tiepoint  coefficients.  Returned  from  USETIE 

c  and  used  in  GETICE.  Refer  to  NODS  Handbook  Appendix 

c  CC  pp  23-24. 

c 

c     Outputs 

c 

c      INTERNAL 

c 

c         outfil  :  The  name  and  descriptor  string  which  reference 

c  the  data  file  which  will  be  used  to  store  the 

c  ice  concentration  grid.  This  string  is  computed 

c  internally  from  the  variables,  indate  and  grdtyp.  The 

c  file  will  be  created  in  the  default  directory. 

c  Program  will  abort  if  file  currently  exists! 

c 

c  The  file  dimensions  will  be  304  1-byte  cell  values 

c  with  448  grid  lines  for  the  grid  type  'N3C  and 

c  316  1-byte  cell  values  with  332  grid  lines  for  the 

c  grid  type  'S3C. 

c 

c        ice      :  Output  buffer  containing  the  computed  ice 

c  concentrations  for  the  25  km  total  ice  concentration 

c  grid.  Totice  is  316  bytes  in  length  which  corresponds 

c  to  the  length  of  an  S3C  raster  line. 

c 

c  Vince  Troisi 

c  5/3/1989 

c 

c  1/4/1989  -  Modified  for  SSM/I  CD-ROMs 

c  7/1/91  -  Modified  for  Apple  Macintosh  and  25  km  grid  by  W.  A.  Wright 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

integer  idiml,  idim2,  idim3 
c 

c  There  are  five  channels  in  the  N3B  or  S3B  grid  cells 
c 

parameter  (idiml  =  5) 

parameter  (maxlen  =  316) 

integer  iunit,  ounit,  munit,  olen,z 
integer  ymd(6),  ier,  recln,ok 

character  *  3  grdtyp,  icetyp 
character  *  9  drive 
character  *  2  year,  month,  day 
character  *  6  indate(95),  ifile 
character  *  3  mon(12) 
character  *80  outfil,infile,  mfile 
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integer    *2  buff(maxlen*idiml),abyte,bbyte 

integer    *  1  ice(maxlen) ,  mask(maxlen),atemp(2),btemp(2) 

real  c(12) 

equivalence  (bbyte,btemp) 

equivalence  (atemp,abyte) 

data    mon 

+     /'JAN'/FEB'/MAR'/APR'/MAY'/JUN'/JUL'/AUG'/SEP', 
+     'OCTVNOVVDEC/ 


C= 


c  Query  user  for  search  criteria 

c  Query  for  name  of  CD-Rom  drive 
print  *,  'Name  of  CD-Rom?' 
read  (*,  '(a)')  drive 

c  Query  user  for  file  dates 

5        print  *,  'Number  of  files  to  create?  (95  max)' 

print  *,  'Note:  Read  lasts  approx.  10  min.  per  file' 
read  *,  n 

c      print  *,  'Enter  file  dates  (yymmdd).' 

open  (unit=29,  file='Archs  HardDisk:THESIS:daylist', 
+         status='old') 
do  10i=l,n 
read  (29,'(a)')  indate(i) 

10       continue 

print  *,  'These  are  the  selected  dates:' 
do20i=l,n 
print  *,  indate(i) 
20    continue 

print  *,  'OK?  (1-yes,  2-no,  3-stop)  <ret>' 
read  *,  ok 

if  (ok.eq.2)  go  to  5 
if  (ok.eq.3)  go  to  6000 

c  Query  user  for  appropriate  polar  region 

print  *,'  ' 

print  *,'Select  one  of  the  polar  regions:' 
print  *,'  1.  North' 
print  *,'  2.  South' 
print  *, 'Enter  1  or  2:' 
read  *,  ihem 
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c  Query  user  for  ice  type  to  compute 

print  *,'  ' 

print  *, 'Select  type  of  ice  characteristic  to  derive:' 
print  *,'  1.  Total  Ice  Concentration' 
print  *,'  2.  Multi-Year  Ice  Concentration' 
print  *, 'Enter  1  or  2:' 
read  *,  icepar 

c  Combine  responses  to  form  filenames  and  paths  and  to  define 
c  10  constants 

if  (ihem.eq.l)  then 

grdtyp  =  'N3B' 

icetyp  =  'N3C 

lrec    =3040 

idim2  =  304 

idim3  =  448 

npar  =  idiml*idim2 
else 

grdtyp  =  'S3B' 

icetyp  =  'S3C 

lrec    =3160 

idim2  =  316 

idim3  =  332 

npar  =  idiml*idim2 
endif 


do5000z=l,n 
ifile  =  indate(z) 
read  (ifile(3:4),'(i2)')mm 

infile=drive//':*//grdtyp//':19'//ifile(l:2)//*:'//mon(mm)// 
+         ,:7/ifile(l:6)//*.'//grdtyp//';l' 

print  *,  'infile  =  ',  infile 

c  Open  brightness  temperature  file  on  CD 

iunit  =  21 

open  (unit=iunit,file=infile,status='old',access='sequential', 
+       form='unformatted',recl=lrec,iostat=ier) 
if  (ier.ne.0)  then 

print  *, infile, '  not  accessible!  <ret>' 
pause 
stop 
endif 
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c  Construct  output  file  name  and  open  the  file. 

outfil  =  'Archs  HardDisk:THESIS:Conc  Data:'//ifile(l:6)// 

+  '.'//icetyp 

print  *,  'output  file  =  ',outfil 
ounit=  22 

open  (unit=ounit,file=outfil,status='new',recl=idim2, 
+    iostat=ier,form='unformatted',  access='sequential') 
if  (ier.ne.O)  then 

print  *, 'Error  opening  file  \outfil(25:35) 
stop 
endif 

c  Open  Landmask  file 

mfile  =  drive//':TOOLS:*//grdtyp//'MASK.DAT;l' 

print  *,  'mask  file  =  ', mfile 

munit  =  18 

open  (unit=munit,file=mfile,  access='sequential', 

+  fonn='unformatted',recl=idim2,status='old') 

c  Get  ice  algorithm  coefficients  derived  from  tiepoints 

call  usetie  (c,ihem) 


c start  ice  concentration  algorithm 

dol0001p=l,idim3 

c  Swap  bytes  in  buffer  file  for  Apple  Macintosh 

do  800inb=l,npar 
read(iunit)  bbyte 

atemp(2)=btemp(  1 ) 
atemp(  1  )=btemp(2) 

buff(inb)=abyte 

800  continue 

read(munit)(mask(inm),inm=  1  ,idim2) 
c  Average  TBs  into  25  km  cells  and  determine  ice  concentrations 

call  ice25  (buff,c,maxlen,idiml,idim2,ice,mask,icepar) 
c  Write  grid  line  of  ice  concentrations  to  file. 

write(ounit)(ice(iout),iout=l,idim2) 
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1000  continue 

c- end  ice  algorithm 

close  (unit=munit) 
close  (unit=iunit) 
close  (unit=ounit) 
close  (unit=29) 
5000  continue 

6000  end 
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program  locate 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

c  LOCATE  -  This  program  transforms  I, J  coordinates  of  an  SSM/I  grid  cell 

c  to  latitude  and  longitude  coordinates.  This  program  provides 

c  the  inverse  functions  as  well.  LOCATE  interfaces  to  the  revised 

c  forms  of  the  subroutines,  MAPXY  and  MAPLL. 
c 

c  User-defined  Paramters: 
c 

c  gtype    :  Integer  supplied  by  the  user  to  describe  one  of  the  three 

c  grid  cell  dimensions  (12.5  km,  25.0  km,  or  50.0  km), 
c 

c  ihem     :  Integer  supplied  by  the  user  to  describe  one  of  the  two 

c  polar  regions  (l=North  ,  2=South) 
c 

c  itrans  :  Integer  supplied  by  the  user  to  describe  the  type  of 

c  transformation  LOCATE  will  perform  (l=I,J-to-Lat,Lon; 

c  2=Lat,Lon-to-I,J) 
c 

c  i,j      :  Integers  supplied  by  the  user  when  itrans  =  1.  These 

c  integers  describe  the  position  of  a  cell  in  an  SSM/I  grid, 
c 

c  lat,lon  :  Reals  supplied  by  the  user  when  itrans  =  2.  These 

c  integers  describe  the  latitude  and  longitude  in  an  SSM/I 

c  grid  which  LOCATE  will  transform  to  an  I,J  grid  cell  position, 

c  Note:  All  latitudes  and  longitudes  must  be  entered  as 

c  positive  numbers! 
c 

c  Internal: 
c 

c  x,y      :  Distance  in  kilometers  from  the  origin  of  the  grid 

c  (ie.,  pole), 
c 

c  alat, 

c  alon     :  Computed  latitude  and  longitude  returned  from  MAPXY. 
c 

c  SGN      :  Sign  of  the  latitude  (positive  =  north  latitude, 

c  negative  =  south  latitude) 
c 

c  delta    :  Meridian  offset  for  the  SSM/I  grids  (0  degrees  for 

c  the  South  Polar  grids;  45  degrees  for  the  North  Polar 

c  grids, 
c 

c  kk       :  Integer  variable  used  for  reorientation  of  the  grid.  The 

c  grid  is  'flipped'  in  the  Y  direction  for  transformations, 
c 

c  SLAT     :  Standard  latitude  for  the  SSM/I  grids  is  70  degrees, 
c 

c  numy     :  Number  of  lines  in  an  SSM/I  grid.  This  attribute  varies 

c  for  each  of  the  six  grids. 
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c 

c      cell     :  Size  of  the  SSM/I  grid  (  12.5  km,  25.0  km,  50.0  km) 

c 

c      xydist  :  Distance  from  the  origin  of  the  grid  in  the  cartesian  plane. 

c  The  x-y  coordinates  for  the  edge  of  the  lower  left  pixel 

c  is  (3850.0,  5350.0)  for  the  northern  grids  and 

c  (3950.0,  3950.0)  for  the  southern  grids. 

c 

c      RE       :  Radius  of  the  earth  in  kilometers. 

c 

c     E       :  Eccentricity  of  the  Hughes  ellipsoid 

c 

c     E2      :  Eccentricity  squared 

c 

c      PI       :  Pi 

c  Written  by  V.J.Troisi  -  January,  1990 

c  Modified  by  W. A.  Wright  -  July,  1 99 1 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


real  SLAT,E,RE,PI 

real  alat,alon,x,y 

integer  ihem 

real  lat,  Ion 

real  SGN,  delta 

integer  numy(2,3) 

real  cell(3),  xydist(2,2) 

data  numy  /  896,  664,  448,  332,  224,  166  / 

data  cell      /  12.5  ,  25.0  ,  50.0/ 

data  xydist    /  3850.0  ,  5350.0  ,  3950.0  ,  3950  / 

SLAT  =  70. 

RE  =  6378.273 

E2  =  .006693883 

PI  =  3.141592654 

E=  sqrt(E2) 
c 

c  Query  for  the  SSM/I  grid  cell  size, 
c 

print  *, 'Enter  the  grid  cell  dimension:' 

print  *,*  1.  12.5  Km' 

print  *,'  2.  25.0  Km' 

print  *,'  3.  50.0  Km' 

read  *,gtype 
c 

c  Query  for  polar  region  of  interest, 
c 

print  *, 'Enter  the  hemisphere  of  interest:' 

print  *,'  1.  North' 

print  *,'  2.  South' 

read  *,ihem 
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c 

c  Define  the  sign  and  meridian  offset  (delta)  for  the  SSM/I  grids. 

c 

if  (ihem.eq.  1)  then 
SGN    =1.0 
delta  =  45. 
else 

SGN    =-1.0 
delta  =  0.0 
endif 
c 

c  Query  for  translation  type, 
c 

print  *, 'Enter  one  of  the  following  transform  functions:' 
print  *,'  1.  Convert  I,J  to  Latitude,  Longitude' 
print  *,'  2.  Convert  Latitude,  Longitude  to  I,J' 
read  *,  itrans 
c 

c  Start  translation 
c 

if  (itrans.eq.l)  then 
c 

c  Obtain  the  I, J  position  of  the  grid  cell  to  transform  to  Latitude 
c  and  Longitude 
c 

print  *, 'Enter  i,j:' 
read  *,i,j 
c 

c  Convert  I, J  pairs  to  x  and  y  distances  from  origin.  The  grid  will  be 
c  'flipped'  in  the  'Y'  direction, 
c 

x=((i- 1  )*cell(gtype))-(xydist(  1  ,ihem)-cell(gtype)/2.) 
kk=numy  (ihem,gtype)-(j- 1 ) 

y=((kk-l)*cell(gtype))-(xydist(2,ihem)-cell(gtype)/2.) 
c 

c  Transform  x  and  y  distances  to  latitude  and  longitude 
c 

call  mapxy  (x,y,alat,alon,SLAT,SGN,E,RE) 
c 

c  Transform  radians  to  degrees, 
c 

alon=alon*180./PI 
alat=alat*180./PI 
alon=alon-delta 
c 

c  Convert  longitude  to  positive  degrees 
c 

if  (alon.le.0.0)  alon=alon+360. 
if  (alon.ge. 360.0)  alon=alon-360. 
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c  Print  the  latitude  and  longitude  for  the  center  of  the  I,J  cell, 
c 

print  *,alat,alon 

pause 
else 
c 

c  Obtain  the  latitude  and  longitude  pair  and  transform  to  cell  where 
c  that  pair  is  located, 
c 

print  *, Enter  latitude  and  longitude  (positive  values):' 

read  *,lat,lon 
c 

c  Transform  degrees  to  radians 
c 

alat=abs(lat)*PI/180. 

alon=(lon+delta)*PI/180. 
c 

c  Transform  latitude  and  longitude  to  x  and  y  distances  from  origin 
c 

call  mapll  (x,y,alat,alon,SLAT,SGN,E,RE) 
c 

c  Convert  x  and  y  distances  from  origin  to  I,J  pair  (iijj) 
c 

ii=nint((x+xydist(l,ihem)-cell(gtype)/2.)/cell(gtype))+l 

jj=nint((y+xydist(2,ihem)-cell(gtype)/2.)/cell(gtype))+l 
c 

c  Flip  grid  orientation  in  the  'Y'  direction 
c 

kk=numy(ihem,gtype)-(jj- 1) 
c 

c  Print  the  I, J  location  of  the  cell, 
c 

print  *,ii,kk 

pause 

endif 
end 
SUBROUTINE  MAPLL  (X,Y,ALAT,ALONG,SLAT,SGN,E,RE) 

Q£* *********************************************************** ********* 

c$ 
c$ 

C$  DESCRIPTION: 

C$ 

C$  This  subroutine  converts  from  geodetic  latitude  and  longitude  to  Polar 

C$  Stereographic  (X,Y)  coordinates  for  the  polar  regions.  The  equations 

C$  are  from  Snyder,  J.  P.,  1982,  Map  Projections  Used  by  the  U.S. 

C$  Geological  Survey,  Geological  Survey  Bulletin  1532,  U.S.  Government 

C$  Printing  Office.  See  JPL  Technical  Memorandum  3349-85- 101  for  further 

C$  details. 

C$ 

c$ 
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C$  ARGUMENTS: 

C$ 

C$  Variable    Type        I/O    Description 

C$ 

C$  ALAT        REAL*4         I      Geodetic  Latitude  (degrees, +90  to -90) 

C$  ALONG       REAL*4         I      Geodetic  Longitude  (degrees,  0  to  360) 

C$  X  REAL*4        O     Polar  Stereographic  X  Coordinate  (km) 

C$  Y  REAL*4        O     Polar  Stereographic  Y  Coordinate  (km) 

C$ 

C$ 

C$  Written  by  C.  S.  Morris  -  April  29,  1985 

C$  Revised  by  C.  S.  Morris  -  December  11,  1985 

C$ 

C$  Revised  by  V.  J.  Troisi  -  January  1990 

C$  SGN  -  provides  hemisphere  dependency  (+/-  1) 

c^******************************************^ 

REAL*4X,Y,ALAT,ALONG,E,E2,CDR,PI,SLAT,MC 
C5> 

c$ 

C$    DEFINITION  OF  CONSTANTS: 
C$ 

C$     Conversion  constant  from  degrees  to  radians  =  57.2957795 1 . 

CDR=57.29577951 

E2=E*E 
C$     Pi=3.141592654. 

PI=3. 141592654 
C$ 

C      Compute  X  and  Y  in  grid  coordinates. 

IF  (ABS(ALAT).LT.PI/2.)  GOTO  250 

X=0.0 

Y=0.0 

GOTO  999 
250  CONTINUE 

T=TAN(PI/4.-ALAT/2.)/((l.-E*SIN(ALAT))/(l.+E*SIN(ALAT)))**(E/2.) 

IF  (ABS(90.-SLAT).LT.l.E-5)  THEN 

RHO=2.*RE*T/((l.+E)**(l.+E)*(l.-E)**(l.-E))**(E/2.) 

ELSE 

SL=SLAT*PI/1 80. 

TC=TAN(PI/4.-SL/2.)/((l.-E*SIN(SL))/(l.+E*SIN(SL)))**(E/2.) 

MC=COS(SL)/SQRT(1.0-E2*(SIN(SL)**2)) 

RHO=RE*MC*T/TC 

END  IF 

Y=-RHO*SGN*COS(SGN*ALONG) 

X=  RHO*SGN*SIN(SGN*ALONG) 
999  CONTINUE 

END 
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SUBROUTINE  MAPXY  (X,Y,ALAT,ALONG,SLAT,SGN,E,RE) 

c$ 
c$ 

C$  DESCRIPTION: 
C$ 

C$  This  subroutine  converts  from  Polar  Stereographic  (X,Y)  coordinates 

C$  to  geodetic  latitude  and  longitude  for  the  polar  regions.  The  equations 

C$  are  from  Snyder,  J.  P.,  1982,  Map  Projections  Used  by  the  U.S. 

C$  Geological  Survey,  Geological  Survey  Bulletin  1532,  U.S.  Government 

C$  Printing  Office.  See  JPL  Technical  Memorandum  3349-85-101  for  further 

C$  details. 

C$ 

c$ 

C$     ARGUMENTS: 

C$ 

C$     Variable    Type        I/O    Description 

C$ 

C$     X  REAL*4        I     Polar  Stereographic  X  Coordinate  (km) 

C$     Y  REAL*4        I     Polar  Stereographic  Y  Coordinate  (km) 

C$    ALAT       REAL*4        O     Geodetic  Latitude  (degrees, +90  to -90) 

C$    ALONG       REAL*4        O     Geodetic  Longitude  (degrees,  0  to  360) 

C$ 

C$ 

C$  Written  by  C.  S.  Morris  -  April  29,  1985 

C$  Revised  by  C.  S.  Morris  -  December  11,  1985 

C$ 

C$  Revised  by  V.  J.  Troisi  -  January  1990 

C$  SGN  -  provide  hemisphere  dependency  (+/-  1) 

C$ 

REAL*4  X,Y,ALAT,ALONG,E,E2,CDR,PI 

C3> 

C$ 

C$     DEFINITION  OF  CONSTANTS : 

C$ 

C$     Conversion  constant  from  degrees  to  radians  =  57.29577951. 

CDR=57.29577951 

E2=E*E 
C$     Pi=3. 141592654. 

PI=3. 141592654 
C$ 

SL  =  SLAT*PI/180. 
200  RHO=SQRT(X**2+Y**2) 

IF  (RHO.GT.0.1)  GOTO  250 

ALAT=90.*SGN 

ALONG=0.0 

GOTO  999 
250CM=COS(SL)/SQRT(1.0-E2*(SIN(SL)**2)) 
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T=TAN((PI/4.0)-(SL/(2.0)))/((  1 .0-E*SIN(SL))/ 
C(  1 .0+E*SIN(SL)))*  *(E/2.0) 
IF  (ABS(SLAT-90.).LT.l.E-5)  THEN 
T=RHO*SQRT((l.+E)**(l.+E)*(l.-E)**(l.-E))/2./RE 
ELSE 

T=RHO*T/(RE*CM) 
END  IF 

CHI=(PI/2.0)-2.0*ATAN(T) 

ALAT=CHI+((E2/2.0)+(5.0*E2**2.0/24.0)+(E2**3.0/12.0))*SIN(2*CHI)+ 
C((7.0*E2**2.0/48.0)+(29.0*E2**3/240.0))*SIN(4X)*CHI)+ 
C(7.0*E2**3.0/120.0)*SIN(6.0*CHI) 
ALAT=SGN*ALAT 
ALONG=ATAN2(SGN*X,-SGN*Y) 
ALONG=SGN*AL0NG 
999  CONTINUE 
END 
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program  u,  v  calc 

*  This  program  calculates  the  u  and  v  components  of  motion  from  a  list  of 

*  user  specified  geodetic  latitude  and  longitude  postions.  Results  are  output 

*  as  an  ASCII  data  file  in  units  of  cm/s.  Appropriate  lines  in  code  must  be 
modified  to  specify  the  filepath  of  input  and  output  files. 


* 


*  Variables:  fname   -  input  file  name 

*  infile  -  input  file  path 

*  ofile   -  output  file  path 

*  lat,lon  -  latitude  and  longitude  expressed  as  real  numbers 

*  mlat    -  mean  latitude  between  two  positions 

*  u,  v     -  u  and  v  components  of  motion 

*  day      -  day  of  position  (integer) 

*  length  -  length  in  days  of  series 

*  tmpd     -  time  between  postions  (in  days) 

* 

character*  10  fname 

character*70  ofile,  infile 

real  lat(60),  lon(60),  u(60),  v(60),mlat 

integer  day(60),  length,  tmpd 

c  Construct  input  file  path  and  open 

5      print  *,  'Enter  input  file  name' 
read  (*,'(a)')  fname 

print  *,  'Enter  number  of  days  in  series' 
read*,  length 

print*,  'Enter  separation  time  period  (days)' 
read*,  tmpd 

infile='Archs  HardDisk:THESIS:Position  Data:'//fname 
open  (unit=15,  file=infile,status='old',iostat=ier) 
if  (ier.ne.O)  then 

print*,infile,'not  accessible!  <ret>' 

pause 

stop 
endif 

c  Construct  output  file  path  and  open 

ofile='Archs  HardDisk:THESIS:Position  Data:'//fname//'.uv' 
open  (unit=18,  file=ofile,  status='new',  form='formatted', 
+        iostat=ier) 
if  (ier.ne.O)then 
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print*, 'Error  opening  ',ofile,'  <ret>' 
pause 
stop 
endif 

c  Read  in  data 

read  (15,10)  (day(i),  lat(i),  lon(i),  i=l, length) 

c  Compute  u,  v  components  (km/time  period) 

do  100  i=l,  length- 1 

v(i)  =  (lat(i+l)-lat(i))*  111.6114 

mlat  =  (lat(i)  +  lat(i+l))/2. 

call  gcircl(mlat,lon(i),mlat,lon(i+l),u(i)) 
if  (lat(i+l).lt.  lat(i))  u(i)  = -u(i) 
print*,  u(i),  v(i) 

100       continue 
pause 

c  Convert  u,  v  components  to  cm/s 

do  200  i=l,  length- 1 

u(i)  =  u(i)  *  (l./tmpd)*(1000.)*(100.)*(l./24.)*(l./3600.) 
v(i)  =  V(i)  *  (l./tmpd)*(1000.)*(100.)*(l./24.)*(l./3600.) 

200       continue 

c  Write  to  file 

write  (18,20)  V,  V 

write  (18,20)  'cm/s',  'cm/s' 

write  (18,30)  (u(i),  v(i),  i=l,length-l) 

10    format  (i6,3x,f5.2,2x,f'6.2) 

20    format  (\x,  a,',',3x,a) 

30    format  (lx,  f6.2,',',3x,f6.2) 

end 
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************************************************************************ 

subroutine  gcircl(latl,lonl,lat2,lon2,gcircl) 

c  This  program  computes  great  circle  distance. 

c  Inputs  are  latitude/longitude  pairs  in  decimal  degrees. 

c  Output  is  circumferential  distance  between  the  points  in  kilometers. 

implicit  real*8  (a-z) 

real* 8  dpr 

real*4  latl,  lonl,  lat2,  lon2,  gcircl 

integer  again 

dpr  =  57.29577951308232 


c  Compute  difference  in  longitudes 

t  =  (lon2-  lonl) /dpr 

c  Compute  r 

dsinr  =  dsin(t)  *  dcos(lat2/dpr) 
r  =  dasin(dsinr) 

c  Compute  k 

dsink  =  dsin(lat2  /  dpr)  /  dcos(r) 
k  =  dasin(dsink) 

c  Compute  angle  D 

dcosd  =  dcos(r)  *  dcos(k  -  (latl  /  dpr)) 
d  =  dacos(dcosd)  *  dpr 

c  Compute  distance  based  on  1 1 1 .61 14  km/deg 

gcircl  =  d*  111.6114 

return 

end 
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subroutine  usetie  (  c  ,  ihem  ) 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

c  USETIE  -  This  subroutine  computes  coefficients  for  the  ice  algorithm 
c  using  the  NASA-approved  tiepoints. 

c 

c    Inputs 
c 

c  c      :  Real  array  of  coefficients.  Refer  to  NODS  User 

c  handbook  for  description  of  coefficients  A-L. 

c 

c  ihem  :  Integer  denoting  hemisphere  (l=north,  2=south).  Used 

c  in  determining  which  set  of  coefficients  to  use. 

c 
c 

c    Internals 
c 

c  TB...W:  Open  water  tiepoints.  One  for  each  of  the  three  channels 

c  used  in  the  ice  concentration  algorithm  (19V,  19H,  37V) 

c 

c  TB...F:  Firstyear  tiepoints.  One  for  each  of  the  three  channels 

c  used  in  the  ice  concentration  algorithm  (19V,  19H,  37V) 

c 

c  TB...M:  Multiyear  tiepoints.  One  for  each  of  the  three  channels 

c  used  in  the  ice  concentration  algorithm  (19V,  19H,  37V) 

c 

c  DW     :  Difference  between  component  tiepoints  for  open  water, 

c 

c  S  W     :  Sum  of  the  component  tiepoints  for  open  water, 

c 

c  DF     :  Difference  between  component  tiepoints  for  first  year  ice. 

c 

c  SF     :  Sum  of  the  component  tiepoints  for  first  year  ice. 

c 

c  DM     :  Difference  between  component  tiepoints  for  multiyear  ice. 

c 

c  SM     :  Sum  of  the  component  tiepoints  for  multiyear  ice. 

c 
c 

c    Outputs 
c 

c  c      :  Array  of  coefficients  which  have  been  computed  from 

c  tiepoints. 

c 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


real  c(12) 
integer  ihem 
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c  Tiepoints  for  open  water,  first- year  ice  and  multi-year  ice 
c  The  first  element  of  each  tiepoint  array  represents  the  tiepoints 
c  for  the  North  Polar  region;  the  second  element  is  applied  to  the 
c  South  Polar  region. 

c   declare  the  two-element  open  water  tiepoints 

real  TB19VW(2),  TB19HW(2),  TB37VW(2) 
c   declare  the  two-element  first-year  tiepoints 

real  TB19VF(2),  TB19HF(2),  TB37VF(2) 
c   declare  the  two-element  multi-year  tiepoints 

real  TB19VM(2),  TB19HM(2),  TB37VM(2) 

real  DW(2) ,  DF(2) ,  DM(2) 
real  SW(2) ,  SF(2)  ,  SM(2) 

c  First  component  of  each  tiepoint  pair  is  for  Northern  Hemisphere 
c  Second  component  of  each  tiepoint  pair  is  of  Southern  Hemisphere 
c  These  tiepoints  approved  for  use  by  the  NASA  Sea  Ice  Algorithm 
c  Working  Group  (1/90) 


data  TB19VW/ 175.3,  175.3/ 
dataTB19HW/  97.7,  97.7/ 
data  TB19VF/ 254.0,  251.2/ 
data  TB19HF/  236.0,  241.7/ 
data  TB19VM  /  223.2,  223.2/ 
data  TBI 9HM  /  203.9,  203.9/ 
data  TB  37  VW/  199.6,  199.6/ 
data  TB37  VF  /  250.0,  248.3/ 
data  TB37VM/ 186.3,  186.3/ 


DW(1)  =  TB19VW(ihem)  -  TB19HW(ihem) 
DF(1)  =  TB19VF(ihem)  -  TB19HF(ihem) 
DM(1)  =  TB19VM(ihem)  -  TB19HM(ihem) 
SW(1)  =  TB19VW(ihem)  +  TB19HW(ihem) 
SF(1)  =  TB19VF(ihem)  +  TB19HF(ihem) 
SM(1)  =  TB19VM(ihem)  +  TB19HM(ihem) 
DW(2)  =  TB37VW(ihem)  -  TB19VW(ihem) 
DF(2)  =  TB37VF(ihem)  -  TB19VF(ihem) 
DM(2)  =  TB37VM(ihem)  -  TB19VM(ihem) 
SW(2)  =  TB37VW(ihem)  +  TB19VW(ihem) 
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SF(2)  =  TB37VF(ihem)  +  TB19VF(ihem) 
SM(2)  =  TB37VM(ihem)  +  TB  19VM(ihem) 

c  Derive  coefficients  from  tiepoints 

c(l)  =  DM(1)*DW(2)  -  DM(2)*DW(1) 

c(2)  =  DM(2)*SW(1)  -  DW(2)*SM(1) 

c(3)  =  DW(1)*SM(2)  -  DM(1)*SW(2) 

c(4)  =  SM(1)*SW(2)  -  SM(2)*SW(1) 

c(5)  =  DF(1)*(DM(2)  -  DW(2))  +  DW(1)*(DF(2)  -  DM(2)) 

1+      DM(1)*(DW(2)-DF(2)) 

c(6)  =  DF(2)*(SM(1)  -  SW(1))  +  DW(2)*(SF(1)  -  SM(1)) 

1+     DM(2)*(SW(1)-SF(1)) 

c(7)  =  DF(1)*(SW(2)  -  SM(2))  +  DW(1)*(SM(2)  -  SF(2)) 

1+     DM(1)*(SF(2)  -  SW(2)) 

c(8)  =  SF(2)*(SW(1)  -  SM(1))  +  SW(2)*(SM(1)  -  SF(1)) 

1+     SM(2)*(SF(1)-SW(1)) 

c(9)  =  DF(2)*DW(1)  -  DF(1)*DW(2) 

c(10)  =  DW(2)*SF(1)  -  DF(2)*SW(1) 

c(ll)  =  SW(2)*DF(1)  -  DW(1)*SF(2) 

c(12)  =  SF(2)*SW(1)  -  SF(1)*SW(2) 


return 
end 


107 


subroutine  ice25  ( tb,c,maxlen,idiml,idim2,ice,mask,icepar  ) 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  ICE25  -  Compute  ice  concentrations  from  NSIDC  CD-ROM  data  at  25  km 
c         resolution, 
c 

c  Software  is  modified  from  the  Vince  Troisi's  original  50  km  version  by 

c  W.  A.  Wright,  Naval  Postgraduate  School, 

c 

c    Algorithm 
c 

c         The  ice  concentration  algorithm  is  an  adaptation  of  the 
c         algorithm  designed  by  the  NODS  Software  Group  (refer  to  the 
c  NODS  User  Handbook,  Appendix  CC,  pp  22-26).  The  algorithm 

c         uses  the  latest  SSM/I  tiepoints  which  were  determined  by  Don 
c         Cavalieri  and  Koni  Steffen. 
c 

c        A  weather  filter  (GR  >  .05)  is  applied  to  correct  for  false 
c        concentrations  computed  in  regions  of  open  water, 
c 

c         Computed  concentrations  greater  than  120%  and  less  than  -20% 
c  are  considered  suspect  and  the  corresponding  cell  is  flagged 

c  as  missing  (-99). 

c 

c    Inputs 
c 

c  tb     :  buffer  of  TB  cell  values.  This  buffer  includes  all 

c  SSM/I  TB  parameters  excluding  the  85GHz  channels, 

c 

c  idiml  :  The  number  of  parameters  in  a  SSM/I  cell  (ie.,  5) 

c 

c  idim2  :  The  number  of  cells  in  a  SSM/I  Tb  grid  line, 

c 

c  maxlen:  The  maximum  length  of  tb 

c 

c  mask  :  A  line  from  the  landmask  for  25  km  grid.  The  mask 

c  line  corresponds  to  the  ice  grid  line  to  be 

c  computed.  (l=land;  0=water) 

c 

c  icepar:  The  ice  concentration  type  to  compute  (l=total, 

c  2=multiyear) 

c 
c 

c    Internals 
c 

c  c      :  Array  of  coefficients  used  to  compute  ice  concentrations, 

c  Refer  to  NODS  Handbook  Volue  II.  Appendix  CC  pp  23-24. 

c 

c  ioff  :  Offset  pointers  used  to  point  to  Tb  values  used  in 

c  this  ice  algorithm  ( 19GHz  V,  19GHz  H,  37GHz  V). 

c 
c  ipos  :  position  in  tb  array  which  is  to  be  summed 
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c 

c  sum    :  Real  register  used  for  summing  four  adjacent  Tb  values 

c 

c  count :  The  number  of  Tb  values  not  flagged  as  missing.  A  Tb  value 

c  of  zero  (0)  is  a  missing  Tb  value. 

c 

c  fyr    :  Integer  register  which  stores  the  calculated  first-year 

c  ice  concentration 

c 

c  myr    :  Integer  register  which  stores  the  calculated  multi-year 

c  ice  concentration 

c 

c  tot    :  Integer  register  which  stores  the  calculated  total  ice 

c  concentration 

c 

c  fyrb  :  Two  -  byte  buffer  equivalenced  to  fyr  register.  Only  the 

c  first  byte  is  the  significant  part  of  the  register. 

c 

c  myrb  :  Two  -  byte  buffer  equivalenced  to  myr  register.  Only  the 

c  first  byte  is  the  significant  part  of  the  register. 

c 

c  GR     :  Gradient  ratio 

c 

c  PR     :  Polarization  ratio 

c 

c    Outputs 

c 

c  ice    :  Line  of  25  km  ice  concentration  grid.  Values  are 

c  represented  as  a  single  byte  percentage  ranging 

c  from  0  -  100%.  Cells  over  land  flagged  with  -88; 

c  cells  for  which  concentrations  were  suspect  or 

c  not  computed  are  flagged  with  -99. 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

integer    idiml,idim2 

integer    maxlen 

integer*2  tb(maxlen*idiml) 

c  Offset  from  start  of  cell  to  three  parameters  used  in  ice 
c  calculations  for  two  lines  of  data 

integer*2  ioff(3) 

integer*2  fyr,myr,tot 

integer*  1  fyrb(2),myrb(2) 

integer*  1  ice(idim2),mask(idim2),fyrice,myrice 

logical    GOOD 

c  Three  TB  parameters  for  ice  concentrations 
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real  BT19V,BT19H,BT37V 

real  c(12) 

real  GR,PR,PRGR 

real  ANF,ANM,DD 

equivalence  (fyr,  fyrb(  1 )) 
equivalence  (myr,  myrb(l)) 

c  Offset  pointers  for  Brightness  temperature  cells  used  in  algorithm 

data  ioff/0, 1,3/ 

c  Begin  processing  buffer 

ii  =  0 

do  1000  i  =  l,idiml*idim2,idiml 

c  Increment  pointer  for  next  Ice  cell. 

ii=ii+l 

c  Calculate  Ice  Concentration 

GOOD  =  .true. 

c  Set  ice  pixel  to  land(-88)  if  over  land 

if  (mask(ii).eq.l)  then 

ice(ii)=-88 
else 

c  Determine  Ice  Concentration  of  the  ice  pixel 

do  80  1=1,3 

ipos=i+ioff(l) 

tbval  =  (float(tb(ipos)))/10. 

if(l.eq.l)BT19V  =  tbval 
if(l.eq.2)  BT19H  =  tbval 
if(l.eq.3)  BT37V  =  tbval 
80  continue 
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c  Set  ice  cell  to  missing  if  any  TB  parameters  are  missing 

if  (BT19V.EQ.0.0  .or.  BT19H.EQ.0.0  .or. 
*  BT37V.EQ.O.O  )  then 

fyrice      =  -99 
myrice      =  -99 
ice(ii)     =  -99 
else 

c  Compute  open  water  weather  filter 

GR  =  (BT37V  -  BT19V)/(BT37V  +  BT19V) 

c 

c  Set  ice  cell  to  0%  concentration  if  weather  filter  >  .05 
c  Filter  value  of  0.05  may  be  altered  at  a  future  date! 

c 

if  (GR  .gt.  .05)  then 
fyrice  =  0 
myrice  =  0 
ice(ii)  =  0 
else 

c  Apply  noniterative  ice  algorithm 

c  Compute  polarization  ration  PR 

PR    =  (BT19V  -  BT19H)/(BT19V  +  BT19H) 
PRGR  =  PR*GR 

ANF  =  c(l)  +  c(2)*PR  +  c(3)*GR  +  c(4)*PRGR 
ANM  =  c(9)  +  c(10)*PR  +  c(ll)*GR  +  c(12)*PRGR 
DD    =  c(5)  +  c(6)*PR  +  c(7)*GR  +  c(8)*PRGR 

c  Compute  First- Year,  Multi-Year,  and  Total  Ice  Concentrations 

fvr  =  ifix(100.0*(ANF/DD)) 
myr  =  ifix(100.0*(ANM/DD)) 
tot  =  fvr  +  mvr 


c  Filter  out  values  outside  the  range  of  -20  to  1 20 
c  This  filter  may  be  altered  at  a  future  date 

c 

if  (fyr  .It.  -20  .or. 

*  myr  .It.  -20  .or. 

*  tot  .It.  -20  .or. 

*  fyr  .gt.  120  .or. 

*  myr  .gt.  120  .or. 

*  tot  .gt.  120)  then 
fyrice      =  -99 
myrice      =  -99 
ice(ii)     =  -99 
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GOOD       =  .false, 
endif 

c  Place  concentration  values  into  byte  data  representations 

if  (GOOD)  then 

fyrice  =  fyrb(2) 

myrice  =  myrb(2) 

if  (icepar.eq.l)  ice(ii)  =  fyrice  +  myrice 

if  (icepar.eq.2)  ice(ii)  =  myrice 

c — 

c  Set  concentration  to  0  if  computed  concentration  between  -20  and  0 

c  Set  concentration  to  100  if  computed  concentration  between  100  and  120 
c 

if  (ice(ii)  .It.  0)    ice(ii)  =  0 
if  (ice(ii)  .gt.  100)  ice(ii)  =  100 
endif 

GOOD  =.true. 
endif 
endif 
endif 
1000    continue 

c  Return  computed  ice  concentration  array,  ice 

return 
end 
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